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FOREWORD 


The Space and Information Systems Division (S&ID) of North American 
Aviation, Inc. (NAA) under Contract AF 30(602)-3638 with the Rome Air Develop- 
ment Center (RADC) of the United States Air Force agreed to perform a 10-month 
study designed to develop digital computer techniques in two areas of interest 
to the RADC tracking facility. First, a differential vorrection geocentric 
orbit computation program for reducing observed data was to be prepared which 
would operate in a near-optimum manner at the RADC computer center, Secondly, 
a computational logic which could be utilized in the tracking process for 
driving the tracking antennae in an open-loop mode was to be prepared. This 
second program would empley general perturbations theory in the definition of 
the predicted trajectory, (This latter task is reported in SID 65-1203-3). 


This report was prepared as partial documentation of the first task. The 
contents present the program logic and FQ@RTRAN listings for the main body of 
the required program, The remaining portion of the Logic for this program is 
associated with the identification, ordering, and smoothing of the raw data; 
this information is presented in the discussion of the processor in a separate 
document (SID 65-1203-2), This document should be reviewed carefully, since 
the interface between the two portions of the differential corrections orbit 
computation program is a magnetic tape containing the smoothed data in a com- 
patible format, 


The differential corrections program (DCP) is designed in such a manner that 
residuals associated with any one of the six types of data 


l. Range 

2. Range-Rate 

3, Azimuth and Elevation 

4, Range and Range~Rate 

5, Range, Azimuth, and Elevation 

6, Range-Rate, Azimuth, and Elevation 


(defined by differencing the observed data and a computed set evaluated on the 
best estimate of the trajectory) are minimized in the sense of minimum variance, 
When this operation is accomplished (recursively), the nominal trajectory is 
adjusted. Thus, after processing data acquired over some as yet unspecified 
interval of time, the trajectory will be known to an accuracy sufficient to 
allow it to be predicted with confidence, At this time, data can be prepared 
to allow the program developed as the second task in this contract to be acti- 
vated to acquire more data. In a sense, then, the system would be self- 
perpetuating. 


Complete discussions of the program logic, and operation of the DCP will 
be presented in the various sections of this document. However, the most informative 
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discussions for the user will probably be included in the analysis of the sample 
problem (section 3.), This conclusion is based on the fact that complete des- 
criptions of program input, program operation, and program output for a real. 
problem are included. The remaining portion of the document will be a necessary 
part of the library of the programmer who wishes to modify or adapt the program 
to another use and of the individual who wishes to fully comprehend the rationale 
which is mechanized, 


This contract has been managed at NAA S&ID by Mr. Js A. Hill and directed 
by Mr. G. E. Townsend. Mr. Townsend also designed the rationale for the 
program, coded the major portion of the logic, performed the preliminary checks 
of the operation, and prepared this document. Final checkout of the program 
was accomplished with the aid of Mr. C. C. DeBilzan. 


The assistance offered by RADC personnel under the direction of Mr. Gordon 
Negus (Program Manager) is gratefully acknowledged. RADC Project 4519 applies. 


This technical report has been reviewed and is approved. 
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ABSTRACT 


This document presents the formulation, computational logic and coding 
information developed for the purpose of effecting the definition of 
geocentric satellite orbits. The rationale for this process is constructed 
around the recursive minimum variance data filter developed by R.E. Kalman 
and a specially prepared magnetic tape generated in the preprocessor 
(SID 65 1203-2). 


The trajectory portion of the program is formulated in the Encke 
manner and includes perturbing accelerations resulting from the first 
3 harmonics of the Earth's potential function, atmospheric drag, solar 
radiation pressure, and solar and lunar gravitation. These accelerations 
are integrated via an uncorrected Gauss-Jackson routine started with a 
fourth order Runge-Kutta process. 
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LIST OF SUBROUTINES AND THETR FUNCTIONS 


DEFINITION 








This routine mechanizes the subroutines 
and lower order driven routines designed 
to generate the trajectory of a satellite 
vehicle by a differential corrections 


INPUE provides all of the data (except 
the observation date) necessary to define 


REED is designed as a simple special 
purpose read package to provide the 
capability for inputting modified station 
data and estimated satellite data. 


A special purpose routine designed to 
assure that the first data card is the 


INTTAL sets up the solution process by 
defining all of the required vectors in 
the computational soordinate frame of 
1950.0 from the information given in the 
true equation of date frame. | 


This routine is designed to provide the 
user the option of reading either the 
satellites position and velocity vectors 
at the initial epoch or the equivalent 


Data for the lunar and soler ephemerides 
and for the 1962 U. S. Standard atmosphere 
are stored in this routine. These data 
are loaded directly into memory at the time 
the program is loaded (BLOCK is not 


A special purpose routine designed to 
DUMP the entire blank COMMON array when 


SECTION ROUTINE 
2el MAIN 
process. 
2olel. INPUT 
the trajectory. 
Betelel RED 
Qeeelolet SIGRD 
scation identification card. 
2e2ole2 INTTAL 
2201.3 POSVEL 
set of orbital elements. 
2.2.1.4 BLOGK 
called). 
20264.8 DADUMP al 
called. 
2e3el. TRAS 


TRAJ mechanizes the working portion of the 
program and serves as the means whereby 

the reference trajectory is computed and 

the traeking stations of the problem checked. 
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DEFINITION 








A TE 


CONIC generates the position and velocity _ 
on an arbitrary conic section defined by fr, 
Y, as a function of time. It is utilized 


for the purpose of defining the Encke refer- 
ence trajectory and the luni~solar 


This routine is designed to function in 
conjunction with CONIC for the purpose of 
iteratively solving Kepler's equation. 


This routine defines time on the conic 
section as a function of a position 


PARTL is the partial derivative of the 
time variable with respect to the anomaly 
variable. This information is utilized in 


MOTION is the driver routine for evaluating 
the total acceleration experienced by a 
vehicle moving in space relative to the 
Encke reference trajectory in the 
computational coordinate frame (1950.0). 


OBLN computes the acceleration resulting 
from the first three zonal harmonics of 


This routine defines the acceleration 
produced by the tenuous atmosphere on a 


ATMS operates in conjunction with DRAG 
and is designed to compute tne instantaneous 
estimate of the atmospheric density. 


SECTION ROUTINE 
Riv dekek CONIC 
ephemerides. 
2.3-1.1.1 SEARCH 
Ci ovlolee TIME 
anomaly. 
ples ealreless: PARTL 
a Newton iteration by SEARCH. 
2.3.1.2 MOTION 
2edele2ed OBLN 
the Earth's potential function. 
AeBel rer DRAG 
spherical satellite. 
2301.22.21 A'TMS 
2.3.1.2.3 ENCKE 


Since the largest contribution to the 
acceleration is included in the reference 
trajectory and since off nominal motion 
produces a modification to the gravity 
vector, ENCKE is required to define the 
change in the acceleration resulting from 
the change in the central force. 
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SECTION 


Seer 


2.3.1.2.4 


2s sekewcled 


Si sekaend 


2.3.1.2.5.1 


2.3.1.2.6 


2.3.1.3 


2.3.1.3.1 


2.3.1.3.2 


2.3.1.3.3 


2.3.1.364 


ROUTINE 





PRESS 


SPOWER 


PERT 


FQ 


EPHEM 


INGRAT 


HSIZE 


START 


DIFTAB 


INTEG 


DEFINITION 





PRESS computes the acceleration resulting 
from solar pressure on an equivalent 
spherical satellite. 


This function operates in conjunction with 
PRESS. Its purpose is to define the ratio 
of the solar power to the speed of light. 


PERT defines the gravitidnal acceleration 
experienced by the vehicle due to the sun 
and moon. 


PERT is formulated in 3 manner similar to 
ENCKE and thus it employs a series to 
evaluate the off nominal nature of the 
acceleration. FQ is that series. 


The position vectors for the sun and moon 
relative to the earth which are required 
by both PRESS and PERT are computed by 
EPHEM from data stored in BLOCK. 


INGRAT is a driver routine for producing 
the first and second integrals of the 
output from MOTION. 


This routine determines the optimum 
stepsize for a Causs—Jackson intregration 
of the equation3 of motion. 


START is a fourth order Runge-Kutta 
integration routine utilized to establish 
the data necesesry to produce the difference 
taole utilized in the Gauss-—Jackson process. 


This routine differences the acceleration 
vectors established oy START and evaluates 
the first and second sums on the leading 


As as 


diagonal of the difference table. 


INTEG mechanizes the Gauss-—Jackson 
integration formulas for a stepwise 
integration of the equations of motion. 
INTEG also steps the leading diagonal 
in the difference table. 
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SECTION 


2.3.1.4 


2.3.1.4.1 


Cade lelen 


2.3.1.4.3 


2-41 


Qebelol 


2.4.1.2 


Zoeleced 


ROUTINE 


TRAK 


UNIT 


EQINOX 


GHA 


FILTER 


KALMAN 


STMAT 


TRANS 


DEFINITION 


TRAK is the driver program which checks the 
tracking stations of the problem and 
determines if any can observe the satellite 
at that time. TRAK transfers to FILTER 

if observation data are available at this 
time. 


UNIT computes the position vector for the 
tracking station and establishes the up, 
east, north unit vectors at the station. 


This routine defines the coordinate 
transformation relating the observational 
coordinate frame (true equator of date) 
and the computational frame (1950.0). 
This transformation is the result of 
nutation and precession. 


GHA defines the Greenwich Hour Angle as 
a function of universal time and the number 
of days since 1950.0, 


FILTER is the driver routine for a Kalman 
estimate of the correction to position and 
velocity (state vector) resulting from 
processing the observed data. 


This routine computes the minimum 
variance estimate of the state vector and 
the covariance matrix for the estimation 
errors. 


The errors at two points on the perturbed 
trajectory are related by employing an 
averaging process -and two estimates of the 


transition matrix obtained using conic 
representation (TRANS)- 


TRANS provides STMAT the conic 
approximations of the transition matrices. 
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SECTION 


2.4.1.262 


2.4.1.3 


2.42.4 


Pehelad 


Pare lied Ws loeb 


24.1.6 


2.5 
205.1 


26562 


ROUTINE 





INVAO 


MEASUR 


ERROR 


UPSTAT 


KLEMEN 


DATAPE 


DEFINITION 








This routine works in conjunction with 
STMAT in the averaging process. It is 
designed to determine the analytic inverse 
of the output of TRANS. 


MEASUR computes the matrix of partial 
derivatives of the observation vector 
with respect to the state vector. 


This routine computes the weighting matrix 
for the Kalman estimate. Data for Station 
errors and noise in the observations are 
employed. 


UPSTAT is designed to determine whether 
the state vector as computed by KALMAN 

is known to a degree which is satisfactory 
to allow the position and velocity to be 
corrected. 


Since the orbital elements of the 
osculating conic are of interest in the 
data reduction problem, ELEMEN is included 
to provide this information at those 

times when data is processed. 


DATAPE defines the time of the next 
observation, the recording station, the 
type of data and the data itself. This 
information is recorded on magnetic tape 
by the preprocessor. 


Mathematical Subroutines 


MATMPY 


MTINV 


Matrix multiplication of two conformable 
matrices in double precision. 


Matrix inverse of a square matrix in double 
precision. 
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SECTION ROUTINE DEFINITION 
Fe Naik eee nO ETI AE Oe CED oy PEE Se Pe Ee ee NENA Ee Oye Re oer ere ee ele se 
ee ee ee ee 
Zeek CHOOSE This routine functions in conjunction 
with MTINV to determine if the matrix 
is singular. 


Asad ADDMAT Matrix addition 

2.5.4 SUBMAT Matrix subtraction 

Ve Fas) TRANSP Matrix transposition 

2.5.6 CROSS Vecto~ cross product 

Deel DOT Vector scalar product 

2.5.8 AMAG Vector magnitude 

eee SINH Hyperbolic sine 

Zebee COSH Hyperbolic cosine 

2.6.3 ARKTNS Computes the single precision arc tangent 


and assigns the angle to the proper quadrant 


2.6.4 DERAQ The Kroneeker delta 
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1. INTRODUCTION 


1.1 PURPOSE 


FS4-305 is a FORTRAN IV IBM 7094 program which was written and checked 
using the standard North American Aviation monitor system (NAASYS - version 13), 
However, the program can also be operated on CDC equipment and other systems 
possessing FORTRAN capabilities, Consequently, to minimize possible system 
incompatibilities, care has been exercised to assure that only the basic fea- 
tures of the system are utilized,- This approach assures that the program will be 
operable on new generation machines (e.g,, the IBM system 360), 


The primary function of this program is to construct an accurate trajectory 
of a ‘geocentric satellite based on a series of observations (e.g., range, 
range-rate, azimuth and elevation, etc.). This task is performed through a 
differential corrections process by generating a reference trajectory and 
minimizing the observed minus computed residuals in a minimum variance sense. 
(The operation of the filter is explained in general in section 1,2 and in 
detail in section 2.4.) 


The trajectory required in the corrections process is generated by numer- 
ical integration of the accelerations relative to a reference conic (Encke's 
method), Accelerations resulting from the off-nominal nature of the motion, 
the Earth's oblateness (first 3 zonal harmonics), solar-lunar gravitation, 
atmospheric drag, and solar radiation pressure are all included. Tracking 
Station data and the corresponding relative position data are generated simul- 
taneously as a function of universal time for the purpose of defining the 
predicted values of the observations from each station. 

The observation data utilized in this program are provided on a magnetic 
tape by the preliminary processor (see FS4-305A, SID 65~1203-2), These data 
have been operated on in two distinct fashions, First, the raw data from each 
tracking station were smoothed by filtering the observations over a restricted 
interval of time (approximately 20 seconds) to a parabolic curve and the mid- 
point of this interval was selected for processing in the differential cor- 
rections program. This operation assures that the effect of random scatter 
in the data can be reduced and assures reasonably efficient operation of the 
differential corrections process by reducing the amount of raw data to be 
considered, Secondly, the data from the various sites were identified as to 
the station which recorded them and the type of information gathered; then 
the complete array was arranged in a chronological fashion to facilitate 
filtering of the data in a recursive mode. 


The interface of the preliminary processor and this progran is provided 
by subroutine DATAPE (data tape). This routine reads the tape and identifies 
the data for use in differential corrections solution, Thus, complete 
consistency is assured, 
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1.2 Program Concept 


The differential corrections orbit computation program is designed 
around the Kalman data filter developed and discussed in Subroutine 
KALMAN, Operation requires that input data be read to identify the stations 
which recorded the observed data, and estimates of the position and velocity 
(or orbital elements) of the vehicle at an arbitrary epoch. At this time, 
the first data point (time, type of data, recording station, and the 
observation vector) is read into memory and the trajectory from the initial 
epoch to the time of the observation is defined. When this is done, the 
position and velocity relative to the observing station are computed and 
an error signal is generated based on this data and the actual information 
recorded. 


The filtering of this information or prediction of a corrected position 
and velocity is accomplished recursively (i.e., one set of observations at 
a time) by employing a minimum variance formulation (see KALMAN). This filter 
has been shown to approximate the true non-linear process very well. Further, 
the filter provides for weighting the data in a general manner not obtained 
with a simple least squares solution. Its inclusion is thus felt to constitute 
near optimum design of the differential corrections orbit computation program. 


After predicting the best estimate of the position arid velocity at the 


first data point, a second data point is read and the process is repeated. 
This operation is demonstrated in the following skesch. 
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read station 
and satellite 
data and time 

















Output of 
Preliminary Processor 
observation data 
(time, station, type 
and data) read by 
DATAPE 


no more data 


predict trajectory 
from time to time 
of data 














compute the numbers 
corresponding to the 
observed data and 
construct the observed 
minus computed residuals 





estimate the corrected 
position and velocity 
by employing the 
minimum variance Kalman 
filter (recursive form) 


An appreciation of the general nature of this program is now available. 


The immediate problem is, thus, to add detail to the flow diagram and 
establish arationale which will assure efficient and accurate operatisn. 
This detail can best be included by dividing the program into logical blocks 
(driver subroutines) which accomplish the desired functions ,and by discussing 
each block in terms of the intended purpose and formulation. Indeed, this 
division has been performed and is outlined fully in the discussion of the 
three principal drivers of the program, MAIN, TRAJ and FILTER. Therefore, 
attention will turn to the first, MAIN. 


es 
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2.0 PROGRAM DISCUSSION 


The formulation, rationale and computational logic for the routines designed 
to perform the differential correction of geocentric satellite trajectories is 
presented in the sections which follow. These discussions attempt to document 
each routine in a complete manner so that subsequent modification to any func- 
tion presently being performed can be facilitated, 


Careful review of these separate discussions is essential to establish the 
complete computational logic for the program and a working knowledge of its 
Structure. However, if interest is centered in the area of program utilization, 
attention should be directed primarily into the structure of the sample problem 
and the associated input formats. (In the latter case, no detailed information 
of the program logic is essential.) 


The mechanism for communicating between the routines of the program is 
the COMMON region DATA. A map of this region and the definition of the 
variables assigned is presented in’ Appendix 1. This information should be 
reviewed in conjunction with any attempt to revise the program logic. 
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2.1 MAIN Program 


Purpose: To mechanize the subroutines and lower order 
driver routines designed to effect the differ- 
ential correction of a satellite's position and 
velocity vectors based on data recorded cn a 
magnetic tape by the preliminary processor 
routine, PROCES, 


Deck Name: MAIN 
Calling Sequence: Nee 
Input/Output : None 
Subroutines Required: INPUT (loads satellite and station data) 


DATAPE (loads observation data) 
TRAJ (generates trajectory) 


Functions Required: None 
Approximate Deck Length: 63 (octal) 
Description and Program Logic: 


MAIN is the driver routine which mechanizes the complete program for the 
purpose of processing observed data for a single (i.e., one at a time) 
identifiable satellite and differentially correcting the estimate of the 
position and velocity vectors describing its motion. In order to accomplish 
this task efficiently and facilitate the development and checkout process, 
the program has been constructed in a series of major blocks (modules), each 
of which is itself subdivided into a number of subroutines. These modules 
serve to provide the input data for the program, to read into memory the 
first observed data point from the data tape, to generate an accurate tra- 
jectory which can be utilized in the computation of the observed minus computed 
residuals, to check the various tracking stations employed in the problem for 
the relative position and velocity information, and to process the observed 
data (recorded on a magnetic tape provided by PROCES) and differentially correct 
the elements of the trajectory. In its simplest form, then, the program can be 
viewed as the link between these modules, i.e., 
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CALL DATAPE ] 





Check Rest of 
Stations 


~ Date Point 


However, due to the fact that data required as input to TRAK and FILTER are 
available in TRAJ and TRAK respectively and are not utilized in the MAIN 
Program for any other purpose, this logic has been modified slightly. 

The modification consists of placing the call statement for TRAK (and FILTER) 
in TRAJ (and TRAK) and constructing a transfer within TRAJ which will assure 
that control is retained by TRAJ until all date have been processed. In 

this sense, ther, TRAJ is the gross data processor and the flow diagram 


reduces to: 


INPUT ) 
CALL DATAPE ) 


TRAT 
(DATA PROCESSOR) , 







aan sane This routine is docu- 

mented with the Filter 
Group since successive 
calls will always come 
following the reduction 
of a given observation. 
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A carefui review of these modules is essential to a thorough aprreciation 
of the functioning of the pregram. For this reason, complete documentation of 
all. routines employed will be presented in the following sections of this docu- 
ment. The task of reviewing this information, however, has been simplified by 
grouping the subroutines when their function relates them to a larger problem 
to assist in the comprehension. 
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2.2 The INPUT Group 


This group of routines is designed for the purpose of providing the com- 
munications link between the users of the differential corrections program 
and the program itself, The group contains « routine which provides nominal 
information to the program (INPUT), a routine which reads specific information 
regarding the satellite (REED), a routine which loads ephemeris and atmosphere 
information (BLOCK), a routine which effects conversion from orbital elements 
to position and velocity if this form of the data is provided (POSVEL), and a 


routine to set up the solution process by initializing various data cells in 
memory (INITAL). 


The group also includes in a sense the routine (DATAPE) which reads the 


data tape prepared by the preliminary processor, However, since this routine 


is a more integral part of the logic associated with the filtering of the 
data, it is documented in the FILTER group, 


St Tees 
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2.2.1 Subroutine INPUT 


Purpose: 


Deck Name: 


’ 


Calling Sequence: 


Math Common, 
Name | Dimension Argument Definition 
- 15 


Input/Output: 


| FORTRAN 
| 1/0 Name 
0 


CON 


0 RE 
RPOIL, 
COEF I 
COEFH 
COEFD 
GMER'TH 
OMEGA 


Q GMMOON 
GMSUN 

0 AU 

0 CONDAY 


CONV 2 






NIN 
NOUT 


bem 


INPUT provides all of the data(except the observation 
data) necessary to effect the differential corrections 
solution of a given trajectory. Data for the physical 
characteristics of the Earth, Sun, Moon, and vehicle, 
station location and error data, position and velocity 
information for the satellite, and an estimate of the 
covariance matrix for errors in the state vector at the 
initial epoch are provided. 


INPUT 


CALL INPUT 


DATA (1) A subarray of common 
containing constants 
for the program 


Re i CON (1) Equatorial and Polar 
Rp 1 CON (2) radii (Km) for the refer- 
J 1 CON (3) ence ellipsoid, the first 
H ae CON (4) 3 coefficients of the 
D 1 CON (5) potential function 
ab 1 CON (6) (Jefferies notation), the 
Ww a CON (7) gravitational constant 
(Km3/sec2), and the spin 
rate for the Earth, 
Am 1 CON (8) Gravitational constants 
Be 1 CON (9) for the moon and sun 
. (Km3/sec*) 
AU 1 CON (10) The astronomical Unit (Km) 
























CON (12) Conversion from days to 
seconds (86400, ) 
Conversion from degrees 


to radiang 





CON (12) 










CON (13) 


Input and output tape 
CON (14) 


numbers 
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FORTRAN Math Common/ 
1/0 Name Name Dimension | Argument Definition 


I SAT - 20 DATA (16) | A subarray of common containing 
data pert*ining to the satellite 
at the initial epoch. 





1/0 | SMASS m I SAT (1) Satellite mass (K,), reference 
AREA A A SAT (2) area (,2), drag coefficient, and 
cD Cp ui SAT (3) surface reflectivity. 
REFLEC R 1 SAT (4) 
: 1/0 | RVEC Yr {or 3 SAT (5) Position and velocity vectors (or 
a,e,i) orbital elements) in the true 
¥ (or 3 SAT (8) equator of date frame at the 
2,w, 6) initial epoch (Km-sec Units). 
T/o | TW (DJ) te 1 SaT (12) The whole and fractional part of 
TF (DJF) 1 SAT (12) the day corresponding to ¥, ¥ 
relative to January 0 1950. 
(1950.0) J.D. 2433282.423. 
I/O | WINDEX - aa SAT (13) Indices utilized to define the 
CHECK ig 1 SAT (14) format of RVEC and VVEC (cartesian 
GO NO be 1 SAT (15) vectors if 1. or orbital elements 
CODUMP a HF . SAT (16) if 2.), whether or not each 
station is to be checked at each 
integration step (yes if 1. or 
no if 2.), whether or not trajec- 
; tory data is to be nrepared for 
each Nth integration step (yes 
if N. or no if 0.), and whether 
COMMON is printed (no if zero). 
0 KCHECK - 1 SAT (17) Fixed point equivalents of SAT 
NOGO é i SAT (18) | (14) and SAT (15). 
I SDA = 250 DATA (36) | A subarray of common containing 
station data. 
: si N (I) ns 10 “ Station identification indices - 


used to select those stations to 
be employed. 





es ae 


once AICS Fw. a a a ee 





FORTRAN | Math Comnon/ 
Name Name Dimension Argument Definition 


STATN SDA (1) Station array (1 to 10) 
; containing position infor- 
K mation for each of the 
Name stations and ea 6 charac- 


ter identification. 


0 HORCOR E 10 SDA (41) Horizon corrections in 
elevation for each of the 
10 stations. 


0 STERR 2 30 SDA (51) Station error array (1 
to 10 stations) containing 
variance information for 
latitude, longitude and 
altitude (_, _, Km@) 


0 SNOISE 2 10 SDA (141) | Station noise array, (1 
to 10 stations) containing 
variance information for 
measurement for errors in 
range, range in rate, 
azimuth on baal 
(Km*, Km@ /sec® ee 


0 Number = 1 SDA (241) | The total number of 
stations employed in the 
tracking program for this 
given simulation 


—- ose — me eo — — — oe ome (eee —e ww — ER TE cee CED ETERS ee ome 


T STT - 105 DATA (286) | A subarray of common 
containing statistical 
information (and related 
data) for the satellite 


1/o | PP P 6X6 STT (70) | Covariance matrix for the 
errors in F and V at to 
(metric units) resulting 
| from smoothing raw data 

as done in the prelimin 


ote td he Lary 
processor. 
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Subroutines Required: RFFD (Relative Read rackave) 
POSVFT (Converts elements to #, % ) 
IVITAT (initializes vectors for prohjen) 


Funetions Required: None 


Approximate Neck 
Tength: 1:52 (octal) 


Description: 


INPUT provides all of the constants, the data for a fixed network 
of ten tracking stations, and the logic for selecting (and modifying) from 
the fixed array of stations those to be employed on the analysis. In 
addition, INPUT serves as the driver for a read routine (REED) which is 
constructed in such a manner that any of the constants previously loaded 
into common can be changed at the time the physical data and estimated 
orbital data for the satellite being tracked are loaded. Upon return from 
REED a test is made to determine whether the 6 components of position and 
velocity or 6 constants of integration in the form of the orbital elements 
(a, e, i,m,w, 0.) were provided. (In the latter case the corresponding 
radius and velocity vectors are computed). INPUT then calls subroutine INITAL 
which initializes the cells in the common array containing vectors used 
in the integration of the "best" trajectory ard in the processing of the 
data. As a final step, all of the input data are printed for reference. 


The first observation which should be made of INPUT is that variance 
p data are given for the station errors (i.e., variance in latitude, longitude, 
and altitude) and for the noise errors (i.e., variance in range, range-rate, 
azimuth, and elevation) under the assurption that the errors are uncorre- 
lated, If calibration of the stations employed in this simulation should 
subsequentiv prove that the rresent formulation is inadequate or could 
be improved, the complete arravs shonld be provided. The incorporation of 
these additional data requires onJy that the dimension of the SM, S®, 
SNOIS®, and STERR arravs be altered and that the FITTER rackage be modified 
accordingly to accept the new data. Allowances for this modification nave 
been made in the cesign of the common array in that extra locations have 
- been provided, 


A map of the COMMON array (DATA) utilized by this routine is presented 
ez Appendix 1 to this repart. This map should be utilized to facilitate 
communication with the program and to aid in anaiysis of the comiunication 
between routines. 


= 15 < 
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The second observation which should be made is that, since the data 
being processed have been smoothed previously to reduce the number of raw 
data points and to remove the random scatter in the raw data, the variances 
utilized in the program should correspond to the smoothed data points rather 
than the raw data. Thus, to prepare these data if raw data variances are 
known, it is necessary to develop the covariance matrix for the difference in 
the true and smoothed values of the observation vector. 


Av= Yorun - Psyoors 
and note that the smoothed values obey the equation 


2 
YsmMoora = L X, Xy “ bd 
; c 
: Mc 


where: Xj, = ty - tytapoint 
a,b,c = coefficients of the parabola which best fits the data. 


Further, the vector: of constants (¢@) which is determined in an unweighted 
least squares fit of the parabola to the data can be represented as 


% = P Vopservea 


(see the preliminary processor documentation of SMOOTH). Now,if the observed 
vector is represented as 


-b ary * 
Yobserved = Yerue + N, 
substitution back into the equation for AV yields 
Pe (I - Mb) Yorye - MON 
= AV onup - KN 


Thus, the desired covariance matrix can be found by taking the expected 
value of the productaY av" . 


BE (AP AP? Yen fA vy Ypl at ~ abgnx® — aiepTat + KAT} 
i) ae 
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But: E (Ym Vat) = 0 
E (Up NT) = 0 
E (N Gt) = 0 
BE (NWT) =v 
Leaving 


B(4¥ aT) = Ky KZ 


which is generally not diagonal even for a diagonal VY, but which (for the types 
of problems attempted during checkout.) can be approximated with a diagonal matrix. 


The operation of INPUT is controlled by data provided from REED and on 
the station identification card (SIC) which must precede all of the card 
‘ input. This is accomplished as follows: the SIC (a card of 10 fixed 
point numbers of FORMAT 10 I 4+) is read and those stations for which an 
index other than zero is recorded are assumed to be utilized in the data 
reduction problem (The numerical sequence of the utilized stations must be 


NN ST 
ree eee ene 


ee 


arrays of numbers provided for these stations and loads these data into 

common for use in subsequent operations (it is thus important to note 

that the remaining station data are discarded and will not be available to 

the program at later time) subject to modification from REED if any or all 
‘ of the data from one or more stations is to be changed. 


Having selected the nominal station data, INPUT calls REED for the 
purpose of obtaining the required satellite data and varifying the nominal 
station arrays. REED provides data to the program by assigning any data 
(previous to the first blank 12 digit field) to the common location 
provided in the first 12 digit field (right justified) of the card. When 
REED recognizes the fixed point number 999 in the first field of a card 
it assumes all numerical data have been provided and begins to look for 
new station names (one to a card) until the location 9999 is read. At 

" this point, all data are assumed to be available and control is returned 
to INPUT. It is important to note that with the exception of the grouping 
of the data cards before 999 and the name cards before 9999, no sequencing 
of the data cards for subroutine REED is required. 





The data expected by INPUT from REED is as follows: 


Date (16) = satellite mass (ke) 

Data (17) = cross section area (m@) 

Data (18) = drag coefficient 

Data (19) = surface reflectivity 

Data (20) = radius vector (X, ¥, Z) (Km) or the orbital elements a, e, 
i in the true equator of date frame (Km, ~ , DEG) 

Data (23) = velocity vector (X, Y, 2) (Km/sec) or the orbital 

‘ elements w,., 9 (true anomaly) (DEG, DEG, DEG) 
Data (26) = whole number of days since 1950.0 to beginning, 0% present 


date 
Data (27) = fractional part of present day in U.T. expressed in days 
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WINDEX - an index utilized to define the option for Data 
(104) ~ (109). WINDEX = 1. if these data are position 
and velocity vectors; and 2. if they are orbital elements 


Data (28) 


Data (29) = CHECK - an index employed in TRAK to bypass checking the 
tracking stations unless tracking data are available at 
the given instant. CHECK = 1. no bypass; = 2. for 
efficiency 

Data (30) = GONO - an index employed in TRAK to bypass checking the 


tracking stations unless the particular station involved 
recorded the data. GONO = 1. no bypass; = 2., bypass 


Data (31) = CODUMP - an index used to assist checkout and as a perma- 
nent record of program constants, etc. If CODUMP = 
nonzero, the entire DATA array is printed following 
initialization. 


Data (356) = covariance matrix for the errors in the initial position, 
and velocity vectors in the coordinate system of rT, and v. 
and in units of Km@ and (KM/sec)@ 


These data can be augmented with station data if desired. This step is 
accomplished by identifying the number of the station to be changed (the 
number of stations on the SIC preceding the station to be changed) and 
loading data into the following locations. 


Data (36 + 4N) 

Data (37 + 4N) 
Data (38 + LN) 
Data (39 + 4N) 
Data (77 + N) 

Data (87 + 3N) 

Data (88 + 3N) 
Data (89 + 3N) 
Data (177 + 4N) 
Data (178 + 4N) 
Data (179 + 4N) 
Data (180 + 4N) 


Latitude in radians 

longitude in radians 

Altitude relative to the reference ellipsoid (KM) 
Station name (6 characters) 
Horizon correction (rad) 
Variance in latitude (rad)© 
Variance in longitude (rad) 
Variance in altitude (KM)@ 
Variance in range (KM)2 
Variance in range-rate (KM/sec)* 
Variance in azimuth (rad) 
Variance in elevation (rad)@ 


ani unrunrea aa 


ee ee 
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Computational Logic: 


load constants 
for Earth, Sun 


and Moon 


Station I.D. Card 






1- - - Nyo station 
indices defining the 
tations to be considered 


READ 





load basic 
data for 


a fixed set 
of 10 stations 








numbers containing 


input elements 








Using Nj- - - Nyo and 
station data, construct 
the compact arrays of 


information for the 
| stations of interest 


> CALL REED }+ 






load satellite data 
and alter the arrays of 
station data if desired 


of orbit a> 
zen 
; ° be Fp, Vp 
Tp Yb 
: establishes position and 
velocity vectors in the 
CALL INITAL frame of 1950.0, and initial- 


Try sVo2E (X, X T) | m,A,CD,R oe 


{ COMMON | 





ie 


izes the vectors required to 
define the conic reference 
trajectory and the state 
transition matrix. 
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2.2.1.1 Subroutine REED 


Purpose: REED is designed as a simple reading routine for numeric L 
data and alphabetical characters, This routine is essen- 
tial to the program if permanent data of the type listed in 
subroutine INPUT are to be retained or changed as desired 

Deck Name: READ 

Calling Sequence: CALL REED 


Input/Output: 




































| FORTRAN | Math Common/ 
ae NIN - aL, DaTA (13) The input tape number. 
I N - 1 - The location in the common 





array (DATA) at which the 
floating point number pre- 
sented in the secona 

field of the card is to he 
placed. This numoer is 
assigned a 12-digit field 
but should be rightly 
justified. 


When N = 999, the outines 
interprets this as a sig- 
nal that no more data 
(numerical) are aveili ble 
and transfers to a second 
mode. In this mode, the 
routine reads one 6 char- 
acter alphabetical, rz-ue 
pear card into memor until 
N = 9999, At this time, 
return is trigperad, 


5 DATA (Ny, Up to 5 consacutive pieces 

N+1, N+2, ef information (floating 

N+3, N+4) point numbers) which are 
te be stored in the common 
array. If any of the 5 
numbers is blank, the bal- 
ance of the card is not 
read. 


- 33 - 
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whether there are alpha-numeric date or not. 


Subroutines Required: None 
Functions Required: None 
Approximate Deck Length: 130 (decimal) 
Description: 


REED is a semi-general purpose read routine which is designed to accept 
cards in an 112, 5E 12.8 Format and assign the data recorded in the second 
through sixth fields of the card to the consecutive locations in common 
(MATA), beginning with the location read in the first field. Several comments 
governing the operation of REED are in order, however, since they affect the 
functioning of REED and thus of the main program. The first is that REED 
will not search a data card past the first blank field which it recognizes. 
(Blanks are interpreted as -0.0; thus, zeros must be input using no sign or 
a plus, e.g., 9.0). If non-consecutive pisces of information are to be placed 
in the program as data or as new constants, it is necessary to start a new 
card following the point of discontinuity (er else fill in the data for the 
intervening location(s)). 


The second observation is that 411 data are assumed to be floating point 
numbers for the sake of simplicity. ‘Thus, if data are to be utilized as 
fixed point numbers, a fixing operation must be performed in the calling 
program after the data have been read. This operation is quite easily 
accomplished by equating a fixed point name to the floating point variable 
and equivalenting both variables to locations in the common array (DATA). 


The third observation pertains to the method utilized to identify the 
time at which data are complete. A signal of quite arbitrary nature can be 


' utilized; however, it was felt that the information should be recorded in 


the first field since this field alone is always read. Further, since the 
DATA array used for common generally requires fewer than 1000 locations, 
it was decided that when N » 999 no more numerical data would be read. A 
card hearing this number muet thus follow all of the data. ~ 
When the program identifies & location WN} 999, return to the calling 
routine is not triggered immediately. Rather, a test is made first to see 
if any alpha-numeric information is to be read. If N< 9999, additional 
cards are read in a format I12, A6. If NY9999, transfer to the calling 
program is effected. Thus, this card must also be present in the deck 


eS ee ee ee oom gene | enumnTare 


——we 


Bh ae 
SID 65-1203-] 


a ni en a mn en meee cee 2 ae “i. ut, Pan ape A 





X = SIGN (X, D(X) | 


IK= KtNel 
DATA(IX) = 0(X) 
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2.2e-leie.l1 Subroutine SICRD 


Purpose: SICRD is a special purpose subroutine which 
operates in conjunction with the program logic 
to assure that the first data card is the 
station I. D. card. If not, execution is 


terminated. 
Deck Name: SICR 
Calling Sequence: CALL SICRD (N) 


Input/Output 


T/O | FORTRAN Math Dimension Conmon/ 
Name Name Argument 
0 N Arg 











Definition 


array of numbers indicating 
which of the ten stations 
of the program are to 

be utilized (1 = yes, 

0 = no) 








input and output tape 
numbers. 


Subroutines Required None 
Functions Required None 
Approximate Deck Length 100 (Octal) 
Discussion: 


This rewtine is designed to read the first data card in an "A~FORMAT" 
and to test the first word on that card to determine if it is in fact the 
station I. D. card (SIC). If the card is identified as the SIC, the words 
in the next 10 fields of the card are tested for zeros. Finally, an arr 
of numbers (x7) is constructed of zéros and ones so that the main program 
will employ the proper network. It is importent to note that since the data 
are read in an A Format any non-zero characters can be employed to identify 
the users desire to include a station. 


If the first card is not the SIC, then execution is terminated with 
an error message. 
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2.2.1.2 Subroutine INITAL 


Purpose: INITAL is designed to set up the solutior to ve ner- 
formed within TRAd by taking the input data and 
storing the position and velocity vectors and the 
associated epoch in cells set aside for the purpose 
of computing the conic reference trajectory and *he 
state transition matrix, 


Deck Name: INIT 


Calling Sequence: CALL INITAL 


Input/Output: 
FORTRAN Math Common/ 
T/O | Name Name | Dimension! Argument Definitior. 
I R, V r, Vv 3,3 SAT (5) Radius and velocity 
SAT (8) vectors in the true equator 
of date frame of refer- 
ence (Km, m/sec) 
1/o | TW, TF be mee SAT (11, 12) | Whole and fractional 
WRK (50, 51) | number of mean solar days 
elapsed since the refer- 
ence epoch of 1950.0 
(JD 2433282.423) 
0 ROTATE NP 3X3 WRK (1) Transformation matrices 
ROTINV pTnT 3X 3 WRK (10) relating true equator of 
date frame of reference 
to the mean equator of 
1950.0 system. 
(rp = NPrso) 
0 RCON ro (0) 3 WRK (28) Position and velocity 
VCON Ve (0) 3 | WRK (31) vectors used to describe 


the conic réefrerence tra- 
Sectory for the Encke 
integration (Km, Km/sec) 










The epoch corresponding 

to the radius and velocity 
vectors used to define 

the conic reference 
trajectory s 







Be A ae 
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Dimension 









etl ee rae ne 


Conmon/ | 
Argument Definition 


WRK (36) Position and velocity 
vectors utilized to define 
the state from which 
errors will propagate 

to the time of data 
acquisition. These 
vectors are in the true 
equator of date frame. 
(Kin, Km/sec) 


The epoch defining the 
point from which errors 
will be propagated 


Position and velocity 
vectors corresponding to 
R and V in the frame of 
1950.0 (Km, Km/sec) 


Position and velocity 
vectors for the satellite 
relative to the conic 
reference trajectory 

(Km, Km/sec) 






Position and velocity 
vectors on the conic 
reference trajectory 

at the epoch t (Xm, Km/sec) 


RCONIC 
VCONIC 






STT (64) The initiel components 
ef the vector composed 
of errors in the radius 
and velocity vectors 
(Km, Km/sec) 
ie a 
Subroutines Required: EQINOX (relates frame of date to 1950.0) 
MATMPY (matrix multiplication) 





Functions Required: None 
Approximate Deck 134 (octal) 
Length: 
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iscussion: 


TNTTAT. is intended to take the position and velocitr vectors in the 
tre equator of date frame of reference and compute the corresponding 
position and velocity vectors in the computational frame (the mean 
equator of 1950.9). This is accomplished by calling FQTNOY and performing 
the following multiplication 


¥e50) = pint Vp 


Having completed this operation, the vectors required to perform the 
jntegration of the trajectory and the computation of the state transition 
matrix are loaded into common, At this point, the process is completed 


by establishing the components of the initial state vector (Lie feEtyyt ) 
v{o . 


a ae 
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CALL MARMPT 
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P5q = T rp 
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2.2.1.3 Subroutine POSVEL 


Purpose: POSVEL accepts data in for form of the orbital elements 
and computes the position and velocity vectors 
corresponding to the epoch of the elements. 


Deck Name: POSVEL 
Calling Sequence: CALL POSVEL (GM A, E, OINC, OMEGAW, OMEGA, THETA, 
R, V) 
auaab. 
FORTRAN | Math Common/ 
Name Name Dimension) Argument Denrere ton 


I CoN (4) The product of Newton's gravita- 
tional constant and the mass of 
the Earth. (Km3/sec2) 


The semi-major axis, eccentricity, 
orbital inclination, argument of 
perigee, longitude of the ascending 
node and true anomaly of epoch for 
the elliptical orbit of interest 

as expressed in the true equator 
and equinox of date frame of 
reference (Km, -, deg, deg, deg, 
deg ) 


0 SAT(S) The position and velocity vectors 
SAT(8) in the true equator and equinox 

of date frame of reference (Km, 

Km/sec) (cartesian coordinates) 





Subroutines Required: None 


Functions Required: SIN, SIND (sine function) 
COS, COSD (cosine function) 


SOR (square root) 
ATAN (are tangent) 
Approximate Deck 231 (decimal) 
Length; 
ee 
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Formulation: 


The orbital elements define a corresponding set of position and 
velocity components (in the same coordinate system). POSVEL is designed 


to compute this equivalent set of numbers by utilizing elliptic fcrmulae. 


However, since these formulae are reasonably well known, they will not be 


= 


developed; rather, those employed will be presented. 


The elements which are assumed to be available make up 2 variation of 
the classic set to be specific. 


a 


e 


Q 


Thus, the 


employing 


To 


Ko 


where ¥y 2 ein™l E . =| 


- 
- 


semi-major axis (Km) 

orbital eccentricity 

orbital inclination to-the true equator of date (deg) 
argument of perigee (deg) 


longitude of the ascending node relative to the true vernal 
equinox of date (deg) 


true anomaly of epoch (deg). Thir element is utilized rather than 
the time of perigee passage due t. he fact that it simplifies 
the resolution process 


dynamics of the problem is introduced through these elements by 
elliptic equations as follows: 


a{i-e®) 


1 +€c0s OQ 


tan-L[e sine, | -90 ¢ ¥, ¢+ 90° 
1 +e cos 5, | 


w 


These quantities fail to fully establish the desired information, however, 
since the orlentstion of the plane of motion has not been introduced, this 
step is accomplisnsed by referring to the following sketches. 


= 52: = 
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> nan 
rerr 
xx" 
> vs a x “a 
V =v sin @ r + v cos s 
“A . 
where T and $ are obtained from 
a ae 
r x x 
aA a ul 
SPE <y = Ty(P)T(1)TA(M) yy 
a ut 
w Zz Zz 
603 PCOSQ-SINGCoSE SIND. COSD SIND. + SING 608 6 COSR. SINC SIND x 
=| “Su CosQ ~cosPCose sin 2 “SING $14. + COSY COS+ COS COS Stare y 
SIL SIMD “SIME COSR case Zz) 


and where T( ) means the transformation constructed by rotating the 
coordinates about the 2 - axis in a cew sense through the 
angle 


T(t) means the transformation constructed by rotating the 
coordinates about the x - axis in a ccw sense through 
the angle i 


Y= 80+ 


ae 
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Compasct i. oneal Logic: 


POSVEL 


eee eee ee, eee eee 


(COS % CosN.-sias yp Cose siasar) x Ra (cos ye suiQ + Susy cose cos.0) ¥ Y + suave sas ge 
C 


3 
$=, aE Se Cos. cose cost sms) X + aed SAID. + COSH Cost ccd ye cosy suse 


ra a(/-e C0 e cos 8) 


v= La(%-“a)]® 


igs tan ‘Le sit 8/( 1+ @ cos 6)/ 


> Aa 
rupep 

wh 

Vevsmxr tvcos 7§ 


RETURN 
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22.1.4 Block DATA 


Purpose: 


Deck Name: 


Calling Sequence: 


Input/Output : 
FORTRAN 

T/O| Name 

0 ALT 1 

ALT 2 

ALT 3 

0 REQT 

% RPOL 
0 STEP 1 
j} STEP 2 

: o | RHOF 

RATE 

0 E PHAM 





BLOCK stores atmospheric density and lapse rate data 
and solar-lunar ephemeris data into cells utilized by 
subroutines ATMS and EPHEM, respectively. BLOCK is 
executed at the time the program is loaded; thus, the 
routine is never called in the program. 


BLOCK 
None 

Math Common/ 

Name |Dimension| Argument Definition 

hy L ATCON (1) the altitudes in Km correspond- 

ho 1 ATCON (2) ing to the ranges over which 

hg 7 ATCON (3) | lapse rate and density data 
for the 1962 U.S. standard 
atmosphere will be quoted. 

R, 1 ATCON (4) the equatorial and polar radii 

Rp 1 ATCON (5) of the earth in Km 

4h, 1 ATCON (6) the intervals at which density 

Ahy 1 ATCON (7) data are quoted (between h, 
and h,, and ho and hz, 
respectively) 

p? 21 TABIE (1) the density and lapse rate for 

K 21 TABIE (22) | the 1962 U.S. standard atmosphere 
at the tabulated intervals 

te 1 EPH (1) the date relative to 1950.0 

stg 1 EPHOM (2) (J.D. 2433282.423) at which the 

At 1 EPHOM (3) tabulated data begin (days), the 

Tm ¥n 270 EPHOM (4) time interval for solar ephemeris 

Ts,Vg! 66 EPHOM (274)|} entries (days), the time interval. 


for lunar ephemeris entries(days) 
position and velocity vectors 
for the moon (Km, Kin/sec), and 


position and velocity vectors 
| for the sun (km, Km/sec) 


| 
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Subroutines Required: None 
Functions Required: None 


Approximste Deck 1 (octal) 
length: 


-~ 60 ~ 
SID 65-~-1203~1 





Le 


daivd 


OLEGAH 1G 
O9CTAO 1G 


DSEING Ws * ee eK KR HR KK KR Ke RK KK KH XK 


94COND 19 


/(e*t=l* (1) WVHd3 ) 


QECONOI * *e ek RK KK RK eK RK KE KR HE KH OK HK KH HE HK HR KK 


OccOrmo 19 
OTEOXO Wd 
Q9eCIW0 719 * 
06zZ0NO 19 * 
oszouole* 
01209 18 * 
6972000 18 
OScuHND 19 
347Z0NO 718° 
o€ZdnO TS 
OZZG%H 18 * 
dizeno td‘ 
09ZONG 19 
J6TIXO 198% 
JSTOXS 79 
ILT9ONO 1G 
I9TONG TE 
OS TOAD 18 
07 TOD 18 
Oc TONO Td 
Oc Tuyo 19 
JTTONO 19 
9OTONO 19 
06V0 40 79 
0800 XO 19 
OLSONO TY 
69009819 
0Ssooxo1d 
O07 IONO TE 
Oc OOD Is 
O0200N0 16 


Ge/scc/Ttt 


stto** 
ULTO* § 
Ssoe j** 
877690° § 


€T-JO7V9°HS 
11T-3S594°T 
O1-38S8°S* 

6-368S°2‘ 


ector* 
S26TO° * 
SeEEec® * 
GIeT°é 


€T-30€%°s* 
T1-3S685°¢‘ 
01-39¢0°8* 

8-39e%°C* 


/ soto** 

oeTa°’ Tvto°* 
sozzo° Z26S20° 5 
86040° * 962S0°* 
Tyst*¢ S2Lt°s 

/ €l-3LES*T* 
CT-3LLS°1* 21-3227 °E § 
T1-3S2676* OT~-38TE*ES 
6-36ST*1* 6-39CB8°1S 
8-F628°6* = L-3¥L6°%/ 


SIY¥S3WSHd3 * + + 8 &® K x 


Viva 


* oo me ek 


eS KR KK RK KH RF OK 


JT1O° 
7STO° 
6c8c0° 
%2690° 


(IZ24THP4 CT) 31LVeIViIVG 


€1T-AS49°? 
CT-1869°9 
OT-JLYVEcy 
6-A76E PS 


(Te t=F*C 1) JOHYI VIVA 


we Re Ke HK Re Ke eH KR K K & * JDUSHDSOWLV * & & & * A RK KR Ke * 


we we KR eK RK KE FX 


x ek F 


we Ke Ke RK KR KR HE KK HE HK HF KR K 


*/99T°BLED/ LOY 


ed3lsS 


*1d3is 


we ow KR KR KR KK KH KK OK 


*Wduy 


* 
* 


/°O0S *°O0T/ 
“7°09O2 


‘1034 


(T2)31Vy 


= ek KX SF 


0°39S6T JO JwWVuYd AHL NI 


ae 02¢ 


©e 


bod 
* 


ed3ais ‘1d3 
**O0T/ € 


we oe eR eK eR He K 


mt NO ot 


ANOS 


LS “/¥8L°9SE9/ Wd 71 


Liv *2idyv 


*Ildy viva 


(6CC)}WVHdS / WOHd3S / NOWWOD 
STL TV /NODLV/ NOWWOD 


itv ‘@1iv 


ee eK 


rc eS 


NOOW 


*€T2)SOHY /318Vi/ 


* 


NOWWOD 


Ke kK Oe 


GNV 


NAS 3HL uOs VIVO ALIDOISJA GNV NOILISOd JOIAOYd GNY 


( 296T - 


CSN] 


QUVONVIS 


s*n } 


AN3SW31L1ViS 33NNDS N43 


3 - 


JYSHdSOWLV GVO SaGuvd 3s3aHl 


ViyG x90719 


430798 
soess 


SID 65-1203-1 


= 61S 


WOooOUO 


YVUUUOYU 


© 


ONUUUOOUOYO 


4 
(" 
i 
\ 
| 





= 


OFL0xO 1d 
OE LOND IG 
OCLONO TS 
OTLONO1S 
002040198 
D690ND 19 
089049139 
OL90NO1E 
9990209 78 
OS9OND1E 
079009 12 
OCIOUND 19 
OZ9ONO Ta 
019040198 
O09ONN IS 
06SONN 19 
O8SONN 19 
OLSOAD IE 
O09SOND1E 
O0SSOHO 18 
Ov7SOnNOTE 
O£SOND 18 
O2SOXN 19 
OTSONO 19 
OOSOND td 
06479ND 19 
Os vONoO WE 
9170019 
097020 39 
CS¥OND 19 
07700 18 
0€VOND 18 
0240019 
OTVOXO 19 
00700 19 
O6t0xO 78 
O8EONO Id 


Se/Ec/TT 


*00-374L0600T °0-*QC-S39SZ0T¥ST°O0-£00 3587224946 °0-*90 ASEHSLZET°O T 


*90 3¥ZSB8TCHE"O £SO JSELLIYCB°O-* TO-382 TE 1096 °O 
*90 3¥2B8BLTE6°O-490 FEZE9SLBT°O *90 JB8E2@772yE°O 


*OO-3TES9EZOT*O T 
“sO 362¢S8SS9€B8°O TI 


/(ESt*vet=l* CLIWVHda } VLVO 


f00-39€EE2SO0E"O 400 ABT9OTS9S°O *O0 3LB8LBLLEL*O- 


#90 38TELTB89IZ°O £90 3A9LZE66BH7°O *O00-FEIBBEOHH "OD 
*00-39Z9€978E°0-*SO 3420222T8°0 $90 34¥7688HT °0 
*O0-3E796Z0L%°O £00 370462L98°O *T0-3870E70E¥°O 
*50 3S6S900SIS°O £90 36ZEveL8E°O *00-3¥8TLEST¥°O 


*9G 31TbZe9¥T°O T 
*00 32S260S08°C 7 
*90 3989909%E°O T 
*90 392ESH9¥7°O-T 
*00 3S8Evt+6LL°S I 


$90 FSESTLSEZS°O *SO 39SC96L68°0~-490 ST69THYST "0-490 306S%662E°O T 


*00-3L798Z26¥2°0 £00 3€9€S820S°0 £00 3689L4TB8°0 
#90 3¥99E7152°0-*90 32860L9T7Z°O *TO-3T8S8l8itc°o 


*90 AOSEZ9VYT*O-T 
*90-3869LCE0T°O T 


*TO 3ZESB84ZOT°O $90 3869ZZ9LT°O-490 Z0E9S7Z27E°9O-*GO ALEOITEIS*S T 


$00-34%I36146722°0-* 00-399LESESE°O-*00 3LO8BLTB86°O 


*90 JALZESOST°O-T 


*9Q 381029462°0-£90 3ECTIOZO¥T* O-400-3914062TH°0-400 372296802 °0-T 
*O0 JZESSZTOL°O £90 ASSHLEZOT°O-*£90 3LZOOSST2Z°0-490 ALETHSE9?*°O-T 


*00-39%60%%84°0-*00 JELELZTN6°0-*00-32ZSSLELTZ°O 


*sO ATS¥IOZOT°O-T 


*GO 3L9¥88BEL°O-*90 FEIOIELGE °0-*00~366L848S7°0-*00 3¥2E96606°0~-T 


*00-346E8ESOE*O-*SO 3ST09"2Z2L°O *SO 3Z20T8O0%L6°O 


*90 3S2Z2T16EE*O-T 


J t€2Tsoo=T* Cl WVHda } Viva 
/ O0O0-3Z69¥SHTE°0-400 A896TLYIL°O-*00 36E88S7272L°0-*90 ALyyZyB2T°O T 
$90 SOBEOZLTZ°O §90 3S999€8S7°0-* 10-375 E882EE °O-*00-JESOZHTZE “O-T 


£00 39S5801096°0-490 3¥96”780LT70 *90 AGSccOsTE°O 


*09-369644S2ZT°O *TO-STEEEZIIE"O £10 AJTLLS6G09NT°O—-*90 3SE6646ST°O 


*90 30ZESO0EE"O *£SO 3S€747S8L1°0 £00-355049S60E°G 
*00 JIB876ESEB°O-*90 39¥LL907T°O £90 320596082°0 


*00-3164L060%°G £00 3E9STOESL*O £00-J3864¥USSH*O-*SO0 309647905 °0 


*90 3692Z6E91°G £90 3A06669ESE*O £00-39€806S24°9 
§10-35696969¥°0-“SO0 306€9Z2962°0-*SO 3ZC9BTESET°O 
*OO-30192SYSE°O £00 36ESETTZ8°O *00~-3SS9D¥FLE°O 
£90 SL86T8YET°O-£90 30¥2S7T9E°O *00-30%T 986°C 


*900 3369S6E€EL°O 490 ZBOTEZOTST°9O-*£90 3L290Z78TL2°9- 


*TO-3L618E862°0 *00-35049ESS2°O ‘00 3688EE5T6°O 


£90 J6L98796T°O-T 


*00-38L626614°0 
*90 32tsezie7e*o 


*00 3€7618718°0 
*90 38SZ2T806¢E°9 
‘gO 389989196°0-T 
*00 AlvI9OTES*O Tt 
*90 3¢69€92S2°0 T 
*90 3€0S18e9T °3-T 


Ast tn 


*90 322L769%7E°O-490 ABGEOEECOT*O £00-3SL HT HZ7ST °O-*£00-399968672T°0-T 
®90 SLLESHE6ES°O $90 FE96H~CSST°O-*90 39LLL179€°0-*SO0 3ZBLLYCOL GT 
/(€9S5=1S (CIIWVHdS 3} ViIVG 


/ 20 300000048°0 ‘£0 3000S2S9€°0 


- (SINSI - IANSWSLVLS 239NDS N33 


: ZL0°6L95 T 


HOBIE 
SOESH 


SID 65-1203-1 


62 


eS 


OlTTHO TY 
DOTTHOIG 
D6ITHAO TVS 
9890THO19 
OLOL¥EI9 
O9UTHO TW 
0S OTN 19 
OvOTHON9 
OcOTHO 19 
O2OTHO 18 
OTOTHD 14 
OST O19 
S6604O 198 
o860NN Ta 
0160019 
09690 TY 
0S603%0 78 
076910 190 
S€690ND 19 
O0260NH 18 
OT60H0 19 
0060xN T¢ 
0680019 
Ossoxo ld 
OL80X0 19 
O980HO 1d 
0S8ONO 18 
o780OND 19 
oesoNNd 
32 800 19 
OT8OxO 19 
G6IBOXO WS 
06200 TE 
O8 LOND IG 
ILLIAO TS 
99192018 
5S LONNIE 


S8/EZ/TI 


“od ATTLEBZET*O—-480 3ST8696S2°U *ON-JILT6USSE*D? £00 
*T9-3S89SO96L°G £SO0 SETIGOZST*U %SO JB9OFHEEBHTO~89G 
*JO-SITELEEVE"O 4G 36608E7SGL°O f£00-326TOLELY*?D *SC 


*92 392622661 °0-£90 J6O0LELZVE"O $9U-FTL97S6972°9 SOC- 


*00 328484182°O $90 308S21S0T°0-*90 3796980TE°G-*90 
*OO-3E8ceusel*O fTC-3T69SGSLZB°O £00 342I€4SS6°0 £90 


*90 AZEZLESIE*O-*GO 360FE9L99°Q §$T0-38LE5990%°S-*00- 


£00 39029S8T6°OD $90 3E79STOST°O-§99 AZLSIB9VE*O-*90 


JypoLceloc°do tT 
JELEGTILY*“O~T 
ADBxVV HBF OD T 
akeercedec*do 1 
JOsessyv71°o-T 
388T28ECE°O-T 
319EC6OCUT*V-T 


/(Cde*oe2=1'(1)IWVHd3 } Viva 


130-329L06€ 92° 0-00 38R84L46E9°0-£90 J3EZ2NGLSIL°O $96 
*90 3LT0EL~9IZ° 0-490 3S4S988S2 °D-*O9-3LLLISISE*O-' 99 
SJI-FITBEVZCLZE*O *SO 3JOBZZOESB*O-*90 394O7ZRLIT°S--*9C 
*39-396€S7617°5-*59 J6TITEOTH°9-*O090-JBSTOLSTI°I-*S0 
*GO ALSS9EZOE*O £90 ZETEZBE9LE °0-*00-34222926E ° 9-490 
£90 JTE9GES9S*°O-*SO JLOCETBIS*O {90 Z3E¥SLIZ6T°O £90 
*90-32 2866582 °0-*00-3988LTIBI7°G-*90 FCRVTIE9E6°D-*90 


£99 3€1T279967°9O *£9G ACEECHOOLI°O-419-3472L65656 °0~*E9- 


‘TO 3L6STEYOT°O-*90 3Z6TSETST°O £90 366S8882E°S ‘S50 
*JO0-38ZO7ZHVET*PO *IO-ABZTBYLS7°C £00 3166L5€EC6°9-'90 
*9) 399162882°O £90 J9CETHZ9T°S $00-326TLZEVETS 899 
899 36€00GE29°O-*90 A¥¥LTZLIT*O *90 3D9ITSE9T*U 490 
*90-362S€6LS64°O £30 3SSL92L76°D *OCU-S9L84I¥ST°G-8SS5 
*GO 3)T9TVE6GE°O $96 3S9SLLSSE°D £90-372ILLI9GH*D £90 
FQO-JELVOETIE*O £50 3S2OZETEH*O-£90 JFEHE7ZO¥I°O-490 


JSO9BTOLT°I-T 
3122729588 °)-T 
JCCIVTESE *U-T 
JALIbLEev2*o-t 
J9I69EIL YL EV-T 
AIZGCOLETE°O-1 
J9OL6SGLTT°O T 
310866154°9 T 
3E29GLECE*0-T 
JSTLOLBEST’O I 
32S994928°9 T 
ayeveT tees T 
3e¢7VO9LLY°O TI 
AL6T7A IGE 70 

AOGSOVBEE "SO T 


t 
1 


/(eve*vet=1* CT IWVHdS } Wiva 


/JO-FETHOIELE*G £90 ABEVSEBIO°O £90 ABLITYZLi°O §90 


*90 SPTE8ZSSZ°O-*S0 3976EHTEZ°O *O0-3SL8ES7ETPEO *G0- 


*90 2,.1T729S86°0 *90 398TE6O9T*O-*99 AT12069Z2¢t °9-* S90 
®TI-37S1L69TS°0- 400-3591 06E12°9-*00 346799896 °O §9u 
*90 398689SZE°0-'90 3LESZEVOT* O-*00-3S4ES6S492 °0~* 90 
69) AVIEVESEL*O $90 Z3TI9LLSY¥T°9O-490 3069E8I9Sz2°9-*90 
*90-369S92286E°0-*90 386605981 °S-*90-36T6Z8S24°O SG 
699 JSTE6SBEET°O—-*99 ASIOLESSE*O-*00-3¥6ITEZS47°9-*D0 
*T9-3)196872T2°0-*40 39TS64Z6e°0-4S) JS9GLEII?*G *9U 
SO09-SEIZT HIS 4*° 0-490 ACSLESSSL*°G-* GU-3ZT89bT 2° O-*S0 
*9) 3962¥LS9T°O $90 396822ESE°0-400-31290€ 0562 °9- £00 
*3) J69697E9L5O-*90 3ANHBLBSGHT°O £99 3eSETSHB7c°G *9U 


- (¢(S)YNS3I - LN3AW3ZLVLS ADYNOS N33 7 


J6UCO¥ETI*I-T1 
JST6TETC?’O IT 
JT666EE68°9 T 
ASSTH99LT° 9-T 
3IT9792ESS°O-T 
30ST900C2°9-T 
JEESESYER°O-T 
3Z28LOeEL 8 °O-T 
J6B71S96E PUT 
AC9TETVLL°O TI 
32E996T0S°9D-T 
JETSOICH?°G-T 


HIOWW 
So0¢Ss 


SID 65-1203-1 


~ 63 - 


a 


_ 3 ao ee 


062140739 
o87c1yo7E 
OLCTAD 38 
o9cT Hoe 
OSzino1E 
Ov72tno719 
oe cl yo 
O22TH9 18 
OT2T NOTE 
OOzT O18 
oO6TT NOTE 
OS TINO 18 
OLITXO 19 
o9TT HONIG 
OSTTNON IWS 
ovTTHp1¢d 
Oc TINO 19E 
O2Ttx0 1d 


G8/EC/TT 


32L996TT?S°O 
3LYOLEZET*O-*80 
3896ETB6C°O 
JL TEcecic*o 
AZvYSEcet°d-*80 
ZECLETBOL°O 
A9ON9TIZIS°O 
JO8OSScet*O--80 
3B9COT8672°O 


JT2SLEBBY"°O 
JTC6LEBLGSC°O 
3L8672TILS*°9-*60 
380029687°0 
AEI~CIEBSZ*°O 
AS TE8OILS*O~-*60 
3609¥%C064°0 
389T7VLESC°O 
327260925 °0-*60 


3eecetect<e°o 
396L98CCT° 0-480 
32020T862°0 
JDS¥VIOET?*°O 
31SS7E8ccTt*O0-*80 
32€060862°0 
3Sev7Seete °O 


32090806%°0 
38B8TO6SZ2°O 
3841909LS°0-*60 
JTLT9SETEY°O 
396%926S2°O 
JIZETIOSGLS*O-*69 
3CLSTOZ64"°O 


39oSETB6C°O 80 
JZ2tvetztc°o *10 
JETEIBZET °O-*80 
ZOL97T86Z°G *80 
3B6TEyzZtZ2°0 ‘10 
328 TS8ZET°O-*80 
agstets6z°o ‘80 
3IT6LOL2TC*°O TO 
aGe2sezcet°o-*s80 


GN3 
3899STILS*O-T 
JTELTZEOBvY°O T 
3SiTZ608Sc°O TI 
S6CESOILS °O-FE 
3As98T66er°oO T 
3S61878S7e°O T 
3SS08092S°0-T 
39S$¢9S964°9 T 
AS64T68S2°0 T 


/(6EE*HVOCHT*(CIIWYHdS ) Viva 


360201862°0 £80 
3S8626212°O *T0 
364>482ET°O-*80 
3912608672 °O *80 
3G2eeTetc*o *t0 
JT ctvecet°*o-*s8o 
3SE18086Z°D *80 


ANAWSLVIS 3399NDS N33 


J£EZ9809L1S °0-T 
3cOVETT67°O T 
AJGL8LitT6Se°o T 
ATSILOILGS°O-T 
ALLTEITEY°O T 
A868L46SGC°O TF 
312%7092S°0-T 


43078 
SOESH 


SID 65-1203-1 


64 


SID 65-1203-1 
ay ae 


*T9900) SI W199 NI HLONST W330 
u 00000 WVHd3 
€2S00 HLONGTI 19339 NIO1ND WOHd3 420718 NOWWOD 
u 900049 @d3isS 
u 50090 Td31S Re) yI006 Wd u €9090 L034 
te) 20000 €Litv uy TooCco cl tv u 00000 TLIWV 
29990 HLONAT 2SI909 NISITYO NQOLV AIDE NOWWOD 
u $2090 jive u 0000) SOHY 
JdAl NOTLIVI0O1 TDEWAS 3dAL NONTLVI04 TOEWAS AdGAL NOTLV301 TOSWAS 
$3999 HIONS1 0939390 NIOTYO a1Vavi XIO TE NOWNDD 
SJIGVIUVA NOWWOD 
Viva w»x30798 
dVW JOVADLS 490179 


G8 39Vd SB/ece/tt GOeS4 


a al RA A 3 RR LENE AMINE, ES PP A AERA NRO RIO 


— 


<2 Lane RA ad RE 





2.2.1.5 SUBROUTINE DADUMP 
Purpose: DADUMP prints the entire contents of the unlabeled 
COMMON region DATA. A complete description of the 
output data may be found in Appendix (COMMON map). 
Deck Name: DADMP 
Calling Sequence: CALI. Di DUMP 


Input /Catput 


FORTRAN 
Name Dimension Definition 
T/O | DATA Unlabeled COMMON region composed 
of the subarrays CON, SAT, SDA, 
STT, WRK. 
CON (17) Output tape number 
Subroutines Required: None 
Functions Required: Nene 


Approximate Deck Length: 533 (octal) 


Discussion: 


The principal purpose of this routine is to provide a print-out of the 
pertinent conditions at program initiation. After reading the program input 
data and performing any necessary initialization, subroutine DADUMP is 
mechanized from subroutine INPUT by setting the flag CODUMP to a-non~zero 
value (input data location 31). 


Utilization of this option provides a permanent record of all of the 
numbers associated with a particular run (gravitational constants, etc.), 
and a convenient means of checking for system jincompatibilities during 
program checkout. 
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2.3 THE TRAJECTORY GROUP 


The second of the major blocks included within the differential 
corrections program is the trajectory program. This group of routines 
evaluates the acceleration vector, integrates the acceleration to define 
the position and velocity and computes the data relative to the selected 
tracking stations. These functions are discharged by four separate driver 
routines. 


1) CONIC (Generates conic reference trajectory for an Encke 
solution to the equations of motion - called from 
INGRAT) 

2) MOTION (Generates the total acceleration vector for motion 


relative to the reference ellipse) 


3) INGRAT (Integrates the acceleration vector obtained from 
MOTION and defines the position and velocity vectors) 


h) TRAK (Defines the position and velocity vectors relative 
to the tracking stations and computes the predicted 
values of the observed data) 


which are mechanized to function in conjunction with a master routine TRAd. 
This master is discussed on the following pages and the manner in which these 
routines are incornorated into the total computational problem 1s demon- 
strated. 
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26361 Subroutine TRAJ 


Purpose: TRAJ mechanizes the working portion of the program 
and serves as the means whereby the reference 
trajectory is computed and the tracking stations of 
the problem checked. Once entered, control is not 
returned to MAIN until the final. data point is 


processed, 
Deck Name: TRAJEC 
Calling Sequence: CALL TRAJ | 


Input/Output : 


FORTRAN | Math Conmon/ | 
whe Name name inensie Argument | Definition 


TCONW beenie WRK(34) [TW and TF corresponding to 
WRK(35) |the epoch at which the 
reference vllipse was 
[defined (days). 



































I WRK(44,) |The position and velocity 
WRK(4'7) | vectors in the frame of 
{1950.0 (Km, Km/sec.) 
I WRK(50) | Whole & fractional part of 
WRK(51) {the number of days since 
1950.0 for present epoch 
(Days from J.D. 2433283.423) 
I AT 3 WRK(52) | Position and velocity 
Av 3 WRK(55) [vectors defining the motion 
relative to the Encke 
| ellipse (Km, Km/sec.) 
pe WRK(104)| time in seconds since last 
reference trajectory 


rectification 






WRK(117) 
WRK(120 


Position and velocity 
vectors at TW plus TF for 
the conic reference tra- 
jectory (Km, Km/sec.) 


RCONIC 
VCONIC 
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Subroutines Required: EPHEM (ephemeris for Sun and Moon position vectors) 


MOTION (driver for evaluating the total acceleration 
experienced by vehicle) 


INGRAT (numerical integration driver) 


RQINOX (routine for relating the true equator of date 
to the mean equator of 1950.0) 


TRAK (routine designed to compute topocentric data 
for tracking stations) 


Functions Required: None 
Approximate Deck 115 octal 
Length: 

Discussion: 


TRAY does not in itself accomplish any portion of the required solution. 
Rather it is constructed for the purpose of sequencing the operations which 
are performed in generating the trajectory and checking the stations for 
visibility. In this sense, then, TRAJ is a direct extension of the program 
entitled MAIN and may appear unnecessary. However, for the purpose of check- 
out, it was deemed desirable to test the input and output functions of the 
program separately. Further, for the purpose of program appreciation and 
understanding it was deemed preferable to retain the checkout structure. 


The first step taken upon entry into TRAJ is the computation of the 
solar and lunar position vectors (EPHEM) to aid in the computation of the 
gravitational and solar pressure contributions to the acceleration vector, 
This done, MOTION is extered for the purpose of computing the total accelera~ 
tion vecter and INGRAT is called for the purpose of estimating the positior 
and velocity vectors st a time subsequent to the present epoch, (the step 
size for this integration is determined by the time interval between two 
successive data points and the magnitude of fifth difference of the accelera~ 
tion vector as described in the discussion of the integration package). At 
this point, all vectors are known at the new epoch, time has been stepped, 
etc, and the conversion to the true equator of date from the mean equator 
of 1950.0 is defined. _ (For the sake of efficiency, this operation is per- 
formed only at one day intervals of time.) As a final step TRAK is called 
to determine if any of the tracking stations being employed observe the 
satellite and if there is data available at this epoch to process, (in the 
latter case, TRAK calls FILTER). 


A transfer is constructed in TRAJ immediately after return from TRAK back 
to the computation of new solar and lunar position vectors. The cycle then 
repeats. Thus, on the surface there appears to be no means for the program 
to terminate. This conclusion is improper since at the time the last data 
point is processed, execution is terminated from within FILTER, 
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Computational Logic: 


CALL EPHEM 


CALL MOTION 


CALL INGRAT 





Time since 

last computation | 
of Nutation and 
precession > 
greater than 

1 day 


CALL EQINOX 
_ CALL TRAK , 


: Transfer as 
' aS data is available 
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2.3.1.1 Conic Reference Motion Group 


The trajectory for the space vehicle is formulated in the Encke 
Manner to assure numerical precision. That is, the difference in the 
acceleration vectors for a vehicle moving on the true trajectory and an 
imaginary vehicle which is moving along a conic trajectory (selected in 
such 98, manner that o5 and ve are equal for the two trajectories) is 
evaluated at every point. This differential acceleration is then 
integrated to obtain the differential position and velocity vectors and 
the true motion constructed by adding the reference position and 
velocities to the corrections. 


This group of routines is designed to provide the reference r 
and v as functions of the time utilizing a deterministic set of variables 
to assure nuwerical significance in all computations. The group is 
completely self contained save for several math functions which are 
provided elsewhere in the program and for the gravitational constant 
which is acquired from common. 
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Purpose: 


Deck Name: 


Calling Sequence: 


Input/Output : 
FORTRAN 
1/o Name 
I RO 
SO 
I TIME 
T/o | xX 


DTIME 





Subroutine CONIC 


CONIC provides the conic reference trajectory to be 
utilized in conjunction with an Encke integration 
being performed in other portions of the program. 

The variables are those recommended by Dr. 5. Herrick 
and are employed to avoid ambiguities in the solution 
as the eccentricity approaches zero or one and/or 

as the inclination of the reference orbit plane 
approaches zero. 


CON 


CALL CQNIC (RG, SB, TIME, X, DTIME, R, V) 


Math Common/ 
Name | Dimension {Argument Definition 


Tae | 3 Arg the position and velocity 
Ving 3 Arg vectors at an arbitrary 


epoch in cartesian coordi~ 
nates (Km, Km/sec.) 


t 1 Arg The time in seconds from 
the epoch of To: V5 to the 
epoch at which r and V are 
desired. 


X 1 Arg The eccentric anomaly 
variables (E ~ Eo) defining 
the position at the desired 
time, The quantity is 
saved to aid in future iter- 
ative solutions for 

X¥=X( t) 

The incremental time in 

seconds since the previous 

solution yielded the value 
of X being fed back into 

CONIC 
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Subroutines Redauired: 


Functions Required: 


Approximate Deck 
Length: 


The radius and velocity 
vectors at time t (in the 
same ) Coordinate frame as 

To, Vo) expressed in Km and | 


Km/sec . 


The gravitational constant 
for ane yo body 
(Km 3/sec@) 





SEARCH (numerical-analytic search for the position 


variable satisfying the time constraint) 


DoT (vector dot product) 
AMAG (vector magnitude) 
COSH (hyperbolic cosine) 
SINH (hyperbolic sine) 
cos (cosine) 

SIN (sine) 

SQRT (square root) 

420 (octal) 
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Formulation: 


The equations of motion for the central force problem are 
Te uli 
Me (1) 


Thus, the rate of change of the angular momentum vectcr is 


dP\e Pree rx(-Ae\r =0 
rae “m eas) (2) 
This equation states that the angular momentum (thus the plane of motion) 
of the resultant conic is constant and leads directly to the fact that any 
vector in the plane can be constructed by linear combination of any two 
vectors selected. i.e., 


P= fie tev 


(3) 
Further, since ¥% and Vo are constants and since V = d/dt (Tr) 
vafargt 
7 (4) 


This representation of r and ¥ is completely definitive for all conic 
trajectories and exhibits none of the ambiguities or indeterminances 
encountered with the solution employing the classic set of elements 

(a, e, i, W,2 , Mo). For this reason, equations 3 and 4 will be employed 
in this program and attention will be turned to the evaluation of the 
quantities f, g, f and g ‘ 


Consider the following sketch which illustrates motion in the orbit 
plane 
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Now since equations 3 and 5 represent the same vector 


a a eats a wn 3) 


(6) 
f> Th vane “S 
‘ (7) 
where hz fap =/u a(/-e') 


The basic problem at this point now reduces itself to the evaluation 
of the angle @ in terms of variables which can be related to time through 
the dynamics of the motion. This task requires that equation 1 be inte- 
grated and that the solution shown to te conic. Howsver, since this process 
has been accomplished in such a large number of references (e.g., Reference 
1), the solution is assumed to be known and the angle @ is recognized to 
be the difference in the true anomalies at positions corresponding to Fro 
ana Yr. 


Under these observaticns the trigonometric functions sin.@ and cos@ 
may be reduced to functions of the eccentric anomaly (selected due to the 


fact that Kepler's equation will be utilized to obtain f, g, f and g as 
function of time) by substituting the following identities: 


sin@ = a yi-e* wan £ (8) 
r 


es [- iS 
coso = rat ae e/ (9) 
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Thus 
aun) #3 = nd (O -2) = wO9 604 ~ sw woe 


xP 


= 9% a (oe es €)- $é, (ce -e] 


«gh =| (seem, - suce)-e(se-s) 
“a }-e Ji | ox- e(se- se,)] 


(10) 
where: the notation s% =sina 3 cx = cos™® andX =E- Ek 
but SE = sin (Eo + X) 

= sXcEo + cX SE 
so that 
van 3 > en /1-2| sx(I-ece,) +€ SE, (-ex+i)] 
ho . 
= 2 ji-€| rv, SX + K'Ve (I-ex)] 
zo - faa 
(12) 
similarly 
cos # = cos@cos & + sin@ sin&® 


_ oa [ (ce -@) (c4,~e) + (/-2*) SE SEs | . 


_ «* [reece + SESE,\+ 07 /~- SE S£o) 
Spee eos 7 


—e(cE +ce,) | 


as) tae Aas re [1-(sxeg, * CX SEo) sé] 
lar ; 


- e(cxce, ~EXSE, + CE.) 
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N 
aR, 


[o* +e? - e* sxce se, - e* cx sz 


“CCXCE, +E SKSE,- ece,] 
(12) 


which can be grouped and reduced as follows 


coo 43 a |(ex- CCXCE,) + CSE, OX + 
K's 


Sx(-e*CE. SE, + CSE.) + (e*)-(ece,) 
3 a | te cx-(&-%) Cx 
Pie tas LQ. 


+ SK( Be 4s) +()- £)- j= ay 


= a fh re (ext!) tp (ey cx+ S(t iM ) 
2 


rr La {aa 


be ve {I+ Eta) ) 
7: rr he nV, (13) 


Coo 8 a? re (Cx+l)+ (BD) (sex) +8x(t 5. ) 
Fr a AA {ua 


22 
—-hM 
AB 


eae) Getfereeala BB) 
Ma 


a j crew! h- Ye) + SX t favo \ 
/ 
7 


a nr ALA e Ma 


nay } (14) 
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. Finally, Kepler's equation in terms of the time of transit from te to 
r can be obtained as 


[Be (é-%) = (£- esé) - (4g - eS E,) 


= X — €( SE ~ S&) 
“XS (SxCk + CK SE, ~ Sp) 


AK sXxX(-€c£&,) ES ED 1 EX) 


~ 


X -(/-E\ sx —- Vo ( /-cx 
Oa a 


in _ 
= (X-Sx) t= SX - 5 J-2@ (15) 
Ren ae ee ea 
and the quantities f and g determined as follows 
fee ( C00 8 - 
Kr 


nev van 4 ) 
lap 


a} 
I 


te) [0-29] Gel - a) «(a 


SF 


ay fa\* fire KV 
Po ar py) See as 
(a) pCa)” Garey 13 


+ TeVe ( I-cx)} 
(Le 


SRG ccs 
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or pect 2 (I-cx). (16) 
similarly 
= +6 a. Vref Sxe hy U- cx]: 
é Va (/-e") oe i 


= fo! [4 Sx + (I-ex)] 
4A a W773 
which may be simplified by observing the form of Kepler's equation to yield 
J 
=/@ m SX £4) -(x-5x)— 8% 
gfe |e sxe fg (44) -(rsx)—g 99) 
= (4-4) —/a2 ( x-sx 
(4-4) —/ar ( x-sx) (27) 


The coefficients f and g can now be determined by differentiating 
Kepler's equation as follows: 


4 (Fe 4) = gfe e sain) £)-(- esate £,) | 


fa’. E(/- 4 coo £) 
Vv a” 
eee 
or E=X* F # (18) 
-~ 88 - 
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and by returning to equation 3 to note that 
~ d 
Pea) 


= CX X 
& i Ve re 


~ _ Va" 5 
ari 


(19) 
and 3 - ¢ 
- fe (rexy(t/e ) 
sas ec (20) 


Equations 3, 4, 15, 16, 17, 19, and 20 now provide the complete 
description of the conic motion problem in a completely deterministic set 
of variables. These equations are utilized as follows to define r and ¥ 
as functions of t - to: 

1) solve 15 for X (iterative) 
2) solve 16, 17, 19, and 20 for f, g, f and é 


3) solve 3 and4 for randv¥ 


At this point Dr. 5. Herrick's notation in "Universal" variables 
(Reference 2) can be adopted by defining the following quantities 


Y 


S=Va $x 
G a (/-cx) 
Ve a." (x-sx) 
T = Vu (é-40) 


Wwe 


= 89 
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re —s 
FoF iE rg’ V 
—s Pal sf ae) 


Vief 2 POM 


and’ £',. -@%, fi, and &' can now be expressed in the form in which they 
will be coded. 


“A 
il 
a 
Vio 


ee = aS 

vd rr 
as 
| a 


It is noted that at every step in the development of these equations the 
implicit assumption of elliptic motion has been made (see equation 8 ). 
However, since for hyperbolic motion 


ah = negative 
F = 4.5 
F-F, = Y= ik 


only the following modifications need be made to provide for a hyperbolic 
reference orbit capability 


es 


Ses 


= a, (1- oad £) 


a2 - Vas (¥-aunb ¥) 

Thus, this generalization will be included. 

1. Townsend, G. E., "Orbital Mechanics” Chapter 3 in The Orbital Flight 
Handbook, NASA SP-33 Part 1, 1964. 


2. Herrick, S., "Universal Variables" published by the author, University 
of California (L.A.) November 2, 1964. 
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2.3.1.1,1 Subroutine SEARCH 


Purposes 


Deck Name: 
Calling Sequence: 


Input/Output: 


SEARCH is the iterative logic utilized to solve 
Kepler's equation for the position variable X 
(see CONIC) as a function of the time elapsed 
since the epoch of rp, Vo- 


LOOK 


CALL SEARCH (TAU, DTAU, X) 


T/O FORTRAN| Math Dimension | Common/ 
___ Nane Name _ Argument Definition 


I TAU T 


oO 


Xx x 


Subroutines Required: 


Function Required: 





Arg Universal time variable 
defined by 7= VUE 


Arg the change in tau since 
the solution was last 
iterated. (Tau for first 
solution) 


WRK(130) | the array of constants 
describing the nature of 
the motion 


Arg the eccentric (hyper~ 
bulic) anomaly variable 
defining position on 4 
conic section as a function 
of t = ty 


None 


PARTL (partial derivative of Keplers 
equation with respect to X) 


~ % ~ 
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Function Required (continued) TIME 


SQRT 
SIGN 


Approximate Deck Length: 190 


Discussion: 


(Time relative to the epoch 
of Yo) Vo) 

square root) 

function for attaching the 
sign of one variable to 
another) 

(absolute value) 


(decimal) 


SEARCH solves a monotonic transcendental equation of one degree of 


freedon in a two stage iteration process. 


First, a guess of the solution 


is made and the function evaluated at the guess and at the guess plus a 
fixed increment. The error in the functions are then computed and the 
value of guess which produced the smallest error in an absolute sense is 
saved. The process then repeats itself until the desired solution is 
approached. At this time, the step size is reduced and the search contin 
ued in stages until the error in the desired root is small enough to assure 
that a Newton type iteration will converge to the answer. At this time 

a variable step size is computed based on the error in the function and 
the slope of the transcendental equation at this argument. 


SEARCH may, thus, be constructed as a general purpose routine for 
functions of this type with provision for incorporating a means of gener- 
ating the initial guess, the transcendental function being solved, and an 
analytic partial of the function with respect to its argument. This 
approach has been taken in the routine being discussed with minor revision 
for the sake of brevity and with only the specific application in mind 
(that of solving Kepler's equation of the anomaly variable as a function 


of time). 
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\ SEARCH 


=< 





Xe=X+A7/Q% 
Xe=47/10an 


[Xe= 7/Q% 
A alXel meo ago" 
/O 





9. = T- TIME (Xe) 








| g,= 7- TIME [Xe tA SIGN(y,)) | 
Aé.l 

ee Oe 
Bal 712) SS 


Xe = Xe tASIGN (9,) 






Xe = Xe 
0 " 21 4 = A/F 
9, * 7-TIME 
0 = PARTL 
Xe =Xet g,/D 
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2.3-1.1.2 Function TIME 


LULPOses 
section. 
Deck Name: TIME 
Calling Sequence: TIME (X) 
Input/Output: 
/O | FORTRAN | Math Dimension| Comnon/ 
Name Name 
I X Xx 1 
I ELEM l. Bl 
De Z 
a 1 
Go 1 
) TIME 1 
Subroutines Required: None 
Functions Required: cos 
SIN 
COSH 
SINH 
SQRT 
* 
Approximate Deck Length: 150 
i 1O3t = 


hatitieniie nis. -s, cok” | eheareee- aoe le matiaaaiatell ee | eaemeemenaatateemmend sommes~ Serie > Sees. 








Definition 


the eccentric or hyperbolic 
anomaly variable defining 
position relative to that 


of an initial epoch 


array of constants used to 
describe the conic section: 


the normalized time 
variabletT (t - t) 


ee 
sine) 

(hyperbolic cosines) 
sie sine) 
square root) 


(decinal) 
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TIME is the function which computes the difference 
in the epochs at two points on an arbitrary conic 
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2.3.1.1.3 Funetion PARTL 

‘ al Purpose: PARTL provides the derivative of Kepler's equation 
8 with respect to the anomaly variable (X). 

Deck Name: PART 


Calling Sequence: PARTL (xX) 





Input / Output: 
ar roxas | Math | Dimension Common / Definition 
Name Neme ___| Argument | 
Z x Xx 4. Arg the eccentric or hyper- 
bolic anomaly variable 
defining position relative 
; to that of an arbitrary 
initiel epoch 
I | EERE es 1 WPK(130) | array of constants used 
Dy 1 to describe the conic 
& 1 section 
Oo | PARTL | 27 cs - the partial derivative of 
: the normalized time variable 
with respect to X 
Subroutines Required: None 
Functions Required: cos (cosine) 
sin (sine) 
cosh (hyperbolic cosine) 
sinh (hyperbolic sine) 
SQRE (square root) 
Approximate Deck Leugtin 120 (decimal) 
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2634142 MOTION and the Motion Group 


This group consists of those routines designed to evaluate the 
acceleration vector in the mean equation of 1950.0 frame (selected to 
eliminate terms in the eauations of motion resulting from rotating 
coordinate systems). Forces resulting from 


the earth's oblateness 

atmospheric drag 

displacement from the reference conic 
solar radiation pressure 

sOlar~lunaz gravitational forces 


WEWD KE 
es ww Sw 


are &1l1 included and the group is designed to be complete to itself in that 
it contains all routines (save general purpose math routines) necessary 

to evaluate the forces in question (i.e., an atmospheric routine, an 
ephemeris routine and a solar power function). These routines will be 
discussed on the following pages. 
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Subroutine MOTION 


Purpose: MOTION serves as the driver routine for all of the pre- 
viously diseussed roitines in the Motion Group (i.¢., 
it is designed to compute the differential acceleration 
vector which will be integrated to define the trajectory 
as a function of time). 


Deck Name: MOTN 
Calling Sequence: Call "OTION (ISTART, INDEX) 
Input/Output: 
FORTRAN Math | Comnon/ 
1/0 Name Name | Dimension Argument Nefinition 
Q RCO Ten 3 WRK (28) the conic position and 
vco Vo s) WRK (31 velocity vectors to 
ee (32) be utilized in subse- 
quent computations of 
the reference trajectorrv 
if the orbit is rectified 
in FNCKF (Vm, Km/sec) 
@) TWCON t 1 WRK on the time of the last 
TFCCN 2 WRK (35 rectification of the 
conic reference (Julian 
date in days) 
I R z 3 WRK tie radius and velocity 
V v 3 WRK (47 vectors on the true 
trajectory (Km and 
| Km/sec) in the frame 
of 1950.0 
ie 1 | WRK (50) whole and fractional 
J. WRK (51) part of the date (days) 
relative to J.D. 2433282,423 
ar 3 WRK (52) displacement, velocity 
ar 3 WRK (55) and acceleration vectors 
Ar 3 (58) rélative to the refer- 







ence conic in the frame 
of 1950.0 (Km,Km/sec) 
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FORTRAN 


Nome 












Common / 


Argument Definition 

















RCONIC WRK (117) | radius and velocity vectors 
defining the present esti- 
mate of the conic reference 


path (Km, Km/sec) 










VCONIC WRK (120) 





ISTART ARG counter which is equal to 1 
set upon rectifying the 
conic reference to restart 


the numerical integration 


counter used to identify the 
source of the call of motion 
(TRAT or START) 


Diseussier: 


Mo diseussion of MOTTON as an indererdent roitine is ecrsiderel essertial 
except insofar as twe points are cennerned, First, vpor exit frer FYCKT, a 
test is made to determine if FNCYF considered the displacement from the conic 
reference to he excessive. Tf so, the conic is rectifi?, the 24 *ferential 
position and velocity vectors are zerred, the epoch is recorded as the time 
of rectification, and the BNCK® acc2leration is zeroed prior to evaluating 
all of the other accelerations, 


The other point is that since the forces due to draz and solar radiation 
pressure are expressed in newtons and since *he accelerations are heing evaln- 
ated ir km/sec? rather than m/sec*, the resultant ferce vector mist be 
qiwvided by 1,000 times the mass to obtain consister* units. 
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Computational Logie: 









CALL ENCKE 





<r R}O OO 


uf 


Veo 
TWCON = TH 
TFCON = TF 


RETURN 
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2.3.1.2.1 Subroutine OBLN* 


Purpose: To compute the non-central nature of the force 


(expressed in the mean equator of 1950.0 frame) 
exerted on the satellite resulting from the second, 
third and fourth spherical harmonics of the Earth's 
potential function. 


Deck Name: OBIN 


Calling Sequence: Call OBLN (ACCEL) 


Input/Output: 


FORTRAN | Math Common/ 
I/O _ Name Nane Dimension Argument | Definition 


0 


* 


Pee] 


ACCEL F 3 Arg the difference in the 
force per unit mass 
exerted by the model 
of the Earth and the 
central force for the 
same satellite position 


RE Re i CON (1) equatorial radius of , 
the Earth (Km) 


CJ,CH,CD |J,H,D Lich CON (3,4,5) | second, third and 
fourth coefficient 
of Jeffrey's poten- 
tial function 


GM 1 CON (6) gravitational constant 
F for the Earth (Km3/sec”) 


AN NP 3°38 WRK (1) the rotational matrix 
relating the frame of 
1950.0 to the true t 
equator of date 


yy 
wo 


RVEC WRK (44) radius vector in the 


frame of 1950.0 





Note: This routine is a version of a routine by the same name originally 
prepared by JPL (one reference is JPL TR 32-223) 
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Suteoutines Required: None 


Functions Required: SQRT 
Approximate Deck 

Length: 250 (decimal) 
Formulation: 


The potential function (seffrey's notation) for a vehicle of unit mass 
moving in the vicinity of an oblate Earth is 





U= -£ f LR (1- 3.ain®L) + HRe (3 -Siein*L) waer L 
r 3 r* r? 
+ DRe ( 3- 30aime?L + 3S 2m “L.) bah 
357" 
where F = r P= position vector in the true equator of date frame 
LL = geadetic latitude = sin" F;) = vin"(F) 


and the force per unit mass may be evaluated by constructing the negative 
gradient of U 


alee al ea | 24) 
ax oY 22 


oats 
However, the force vector (F) would be expressed in the true equator of date 
frame rather than in the desired frame of 1950.0. Thus, before constructing 
the partials, it is noted that the radius rector in the two frames are related 


as follows: 
ade ~ a, Qe Gis _> 
PSA = Cte, Gag 23 R 
431 932 Gay 


and thus that 
a (% RK t Aso Rat Qs, R; ) 


With these substitutions established, the force vector (F . expressed in 
the frame of 1950,0 is 


Fo=-(2u , 2 ,w 
ae ° 7H, 7 aR, 
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Where 


2 
~ eg ai\(3- 782) B uy +tf-2 4 3Z)) 4 
~ De a 





Computational 
Logic: 






compute radius vector in true equa~ 
tor of date 


compute constants for use in evalu- 
ating components of acceleration 
(e.g., Sin L, sin; , Ro/?'s etc, 









construct functions of L to be 
attatched to each of the coeffi- 
cients of the potential function 
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2.3.1.2.2 Subroutine DRAG 


Purpose: 


Deck Name: 


Calling Sequence: 


Input/Output: 




















FORTRAN 
Name 


ROTATE 


ROTINV 


To evaluate the force acting on a satellite (in the mean 
equator of 1950.0 frame) resulting from passage through 
a rotating oblate atmosphere. 

DRAG 


Call DRAG (DF) 


Math Common/ 
Name | Dimension Argument Definition 


a 1 CON (7) Spin rate of the Earth 
A 1 SAT (2) Reference area 
Cp 1 SAT (3) Drag coefficient 
Tol 3x3 WRK (1) Rotational matrices 
relating vectors in 
T 3.303 WRK (10) the true equator 
of date frame to 1950.0 
tT 3 WRK (44) Radius vector (1950.0) 
v 3 WRK (47) Velocity vector (1950.0) 
t i WRK (50) Day number relative 





| to 1950.0 (used as in- 


pub into atmospheric 
density routine) 


Drag force (if A =scuare| 
meters, density =Kg/_ | 
m3, and velocity = Km/ 
sec, then D = newtons) 


Subroutines Required: ATMS (atmospheric density), CROSS (cross-product) 


Functions Required: 


Approximate Deck 
Length: 


AMAG (vector magnitude) 


1h6 (octal) 
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Formulation: 


The drag force is computed in the frame of 1950,0 as follows 
—_/ te 
@ = [T) @ 


where [T] = rotational transformation relating the true equator of date to 
the mean equator of 1950,0 
a = spin vector of the Earth in the true equator of date frame 
w= spin vector in the frame of 1950.0 
Vig =O XT 
Vp = V- Vin 
where r = vadius vector (1950.0) 
V = velocity vector (1950,0) 
V, = velocity relative to the wind (1950,0) 
Vy = velocity of the wind (1950,0) 
B=-20(8T) CoA Y/Y, 
where D = drag force in newtons 
ey T) = mass density of the atmosphere 


CyA drag coefficient times reference area 
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Computational Logic: 


convert spin vector from frame of 
date to mean equator 1950,0 


compute velocity of wind in frame of 
1950,0 


compute relative velocity vector 
(1950,0) 


call atmospheric density (1962 U.S, 
standard with solar activity varia- 
jtion imposed) 


2 


D = 1/2 CyAP V 





RETURN 


END 
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S 4es 2634162.2.1 Subroutine ATMS 


Purpose: ATMS computes an approximate atmospheric mass density 
(Kg/m3) at any altitude (h) between 100 and 700 Km. 
(Below 100 Km the density is set equal to that of 
100 Km and an error message is recorded; above 700 Km 
the density is set equal to zero). The U. S. Standard 
Atmosphere 1962, is used as a reference for a table 
(stored in memory within the Data Input Group) with 
21 entries (10 Km steps 100 <h< 200, 50 Km steps 
200<h <¢ 700). A dichotomic interpolation is then used 
in conjunction with the assumed exponential atmosphere 
(in the neighborhoods of the tabulated point) to pro- 
duce an interpolated density (2 to 3 place accuracy 
throughout the table). These data are then corrected 
for the ll-year sunspot activity cycle. 


Deck Name ATMS 


Calling Sequence: CALL ATMS (RVEC, TJD, DENS) 





Input/Output: 
| Math Common/ 
Name {Dimension Argument Definition 
yr 3 Arg radiue vector in the 
true equator of date 
frame 
+ iL Arg dummy variable (such as} 
the Julian date) which 
can be used to compute 
corrections to 
0 DENS Pp 1 Arg atmospheric mass 
density (Kg/m3) 
I AIT1 | hy lL ATCON (1)| Altitude limits for 
ALT 2 hy 1 ATCON (2)j| density data 
ALT 3 ha 1 ATCON (3) 
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Common/ 
Dimension | Argument Definition 





REQT ATCON (4) | equatorial and polar 
RPOL ATCON (5) radii for the Earth 
STEP 1 ATCON (6) |altitude increment for 
STEP 2 ATCON (7) jdensity and lapse 
rate tables 
RHOF TABLE (1) |density and lapse rate 
RATE t TABLE (22) |at altitudes between 
hz, and h3 
Subroutines Required: None 
Functions Required: SQRT 
EXP 
SIN 
ALOG 
AMAG 
Deck Length: 00524 (octal) 
Formulation: 


Atmospheric density will be approximated utilizing data taken from the 
1962 U. S. Standard Atmosphere and assumed exponential behavior with a 
Single parameter, altitude. Thus, the first task is to compute the radial 
distance from the center of the Earth to the subsatellite point (the inter- 
section of the radius vector to the satellite with the Earth's surface) 
which will be defined to be r. This task, in tur, will be accomplished 
by employing the oblate spheroid representation of the Earth. 


Let R be the satellite position vector and c be the proj.ction in the 
plane of the true equator of date 
r= r* + ro° + 13” 


oF = rn? + ro” 


If c = 0, the satellite is directly above either the North or South Pole, 
and 


“A — 
r = 


133 
SID 65-1203-1 


where, rp = polar Earth radius and ’ = local Earth radius. However, 
if c #0; 











Berl x44 2 
7 (i) > 
Local 
= yJ/7 GE rp sb meridian 
Le asl 
where tanE = A 3 
Cc 
Then from the equation for an ellipse, 
Me pee 
at ig 3 
& 2 2 
aa + Canm~£ = 4, 
or se B 
a + tan* E 
Substituting into (1) 
A= Ly / + Cam*E 
a t Lon*E 
and the height above the Earth is 
h=r-f 


Now the atmosphere is assumed to be constructed in uniform concentric 
shells around the (oblate spheroid) Earth with densities which agree with 
those predicted by the 1962 U. S. Standard Model. For simplicity, a two 
part table of these data has been prepared. Part one covers altitudes 
from 100 Km to 200 Km in 10 Km steps while part 2 covers from 200 Km to 
700 Km altitude in 50 Km steps (there are, therefore, 21 entries in the 
table). Since this region contains all altitudes for which the tenuous 
atmosphere should be considered in studies of satellite motion, the table 
was not extended. Rather, below 100 Km the density is set equal to the 
value at 100 Km, and above 700 Km the density is set equal to zero. 
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For any altitude, h, between 100 Km and 700 Km the values in the table 
corresponding to the nearest altitude above and below h, must be determined 
and interpolated to find the density at h. Let m be the index corresponding 
to the nearest tabulated altitude below h, then, 


(1) between 100 and 200 Km, m = greatest integer te +1 





(2) between 200 and 700 Km, m = greatest integer (4-g00 e447 


Let Ah be the distance between h and the nearest tabulated altitude below 
h, 


(1) if 100 Km ¢h 4 200 Km 


X 
Ah 


tt i 


greatest integer in dine f 
(h-100) - X (10) 


(2) if 200 Km ¢h ¢ 700 Km 


X 


greatest integer in (2 f200 sd) 
Ah 


(nh-200) - X (50) 


tol 


The tabular values of density above and below h can now be extrapolated 
to predict two values of the density at h; then, an interpolation can be 
made between the two predicted values of density to include the non- 
exponential nature of the atmosphere. 


Let Ah be the distance between h and the nearest tabulated altitude 
above h, then 


Ah = Ah ~- 10 if 100 h € 200 
Oh = Ah - 50 if 2004 h 4 700 


and let», = extrapolation of lower value (index m) and 72 = extrapolation 
of upper value (index m + 1) 


-Am 4h 
= Ao» e ™ 


-(emri g 
2 = Pans & means 


~ 


where I 


fl 


tabulated value of density 


r tabulated slope (lapse rate) = VA 


Finally, by linear interpolation between 
= bent (1-£)A 
= alee gu ae 
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Computational Logic: 












Ts Rs Cc 
F, = Vv/+ 
Fa “(Bey + 7 


Re» Re (F/ra) _| 








= HA, amd 1 <A, = 
Ro LAD 
(RETURAN [Nd > Az and HS Aye 


H=h-A, 


Hs= 


Ci 3 A,’- A: 

a C,/s / 
(4c hs + Ca +/ 
XM: H/Hs 


C2: 
MM? 
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KC = ¥m 
XM= KC 
Ah= h-XM* Hs 

A2 = Ah- Hs 

A} = Mtl 

Pi * Pg (M) Leo (“A(m)AR) 
2 = Pp (N) “9 (9 (ay Ad) 
xh= Ah Hs 

A> Nh Agr (l- xh) e, 
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26301268 Subroutine ENCKE 


Purpose: 


Deck Name: 
Calling Sequence; 


Input/Output: 


ENCKE provides the correction to the acceleration vector 
resulting from motion on the true (rather than the 
reference) trajectory 

ENCK 


CALL ENCKE (ACCEL, IRECT) 


, FORTRAN Math Common/ 
1/0 Name Name |Dimension| Argument | Definition 
I. GM we 


Gravitational constant 
for the central body 





I RCONIC ' | WRK (117) Position vector on the 





conic reference 
trajectory 


Displacement relative 
to conic position 
(AreT -p 
Differential accelera~ 


tion vector 


Signal for MOTION to 
alter tne reference tra- 
jectory 


Subroutine Required: .None 


Functions Required: 


Approximate Deck 
Length; 


DOT ee product ) 
SQRT (square root) 


213 (octal) 
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Formulation: 


Encke's method of integration involves the computation of the accelera- 
tion vector relative to a known reference which in the present case is 
an ellipse. Thus in addition to real perturbing forces there are terms 
which must be included for the off-nominal nature of the motion. Consider 


Peery F (1) 
ra 
aati (2) 
where += radius vector (in some frame) for true trajectory 


cnet 
F = summation of all nonecentral forces 


P = radius vector (at the time corresponding to the 
instantaneous position on the true trajectory) on the 
conic reference 


definin, oe ae? oe 
: Ar=ry (3) 
The differential equations for AF can be obtained 


Ar =F tu) Aa- | 


U4 


But since lar! may be small, the term in the brackets may be known to 

- less accuracy than F or * due to the fact that two nearly equal numbers 
are differenced. Thus, a modification to maintain accuracy for such 
cases will be presented. 


=F * ae - -(1-fR) F] fp! ai 
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Now fq has one more degree of freedom than does the ratio (7, )3 SO the fol- 
lowing identity can be written arbitrarily 


(1efq) = (142q)°3/2 = (fr)? (5) 


. fy t-(tragy* (a) 


= - 3. a me 5° 7: ¥ 
og A gt 8899" BELT G 


R. 3. 
# (6b) 
where 
2 
2/2 
ry eae | 
= 707 
Z £. 2 (7) 
fe 
/ a a —h > > 
= 2 {(4+at) (2 + SY) “A +0 
ie _ 
= yo: Ar ‘RAL Ar 
af 
Equation 6b,when truncated at any given term,defines the accuracy to 
be obtained from equation 4 for the case where /7/7~ /,0 . Thus, 


if a fixed number of terms are consistently carried, the maximum value 
that q can attain before the reference trajectory is altered (rectified) 
{or before the equation for evaluating fg is changed from equation 6b to 
6a) is defined by the maximum error allowable in the scalar fq. For 
example, if this allowable error is arbitrarily set at 10-8, and if 
seven terms in fq are carried, the upper limit in q is defined by 


t 
3+5+7+9+11+13-15 day < 1078 
D+ 3h 5-67 
rule = 2.5 xX loll 


which corresponds roughly to 
Qmax = +03 
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Computational Logic: 


ENCKE 


compute q, (/r)3, fy 






qgq=(4-Ar +43 ar -Ar ) foe 


(P/r)3 = (1 + 2a) 


compute ~y 


Are B81 (%g) 4 -(2)° 47 





g 2:03 
/RECT =2 
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2.3.1.2.4 Subroutine PRESS 


Purpose: To compute the force exerted on a spherical satellite 
due to the pressure of the impinging light 


Deck Name: PRES 


Calling Sequence: Call PRESS (SFOR) 





Input/Output: 
I/o» | FORTRAN Math Common/ 
Name Name Dimension Argument Definition 
0 SFOR 3 Arg Solar pressure force 
vector in the firame 
of 1950.0 
I RVEC 3 WRK (44) radius vector in the 
frame of 1950.0 
i RSUN 3 WRK (98) Position vector for 
the sun relative to 
the earth (1950.0) 
I A ns SAT (2) Cross sectional area 
of spherical satellite 
I R 1 SAT (4) Surface reflectivity 
I RE ui CON (1) Equatorial radius of 
‘| the earth 
Subroutines Required: None 
Functions Required: AMAG (vector magnitude) 


DOT (Dot product) 
SPOWER (solar power) 
SQRT (square root) 


Approximate Deck 
Length: 220 (octal) 
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Formulation: 


The solar force exerted on a spherical satellite resulting from 
impinging light can be obtained directly from the expression relating the 
solar pressure, i.e., 


is the power of the incident radiation 

is the speed of light 

is the surface reflectivity 

is the angle between the surface normal and the 
impinging light. 


where 


RWartd 


This is accomplished by integrating the pressure over the surface and 
(due to symmetry) applying the result in the direction of the impinging 
light. 


F< P (uu) { cos’ & aA 
> C A 


But for the spherical satellite 
dA = 27 S* sinade 


where S =radius of the satellite 
a =angle between an arbitrary radius vector to the 
satellite skin from its center and the vector 
defining the direction of the impinging light. 


and the integration over the surface reduces to the integration with 
respect to % from 0 to 7/2. The result is 


Fg =z P (14R) 2 S® 
C 3 
Fy= -Fy (fs - F) /|Ts - 7 


This description of the problem is not, however, complete since it 
fails to consider the effects of any time which might be spent in the 
earths shadow. This correction might be made by zeroing P for these 
times but for reasons of data availability in PRESS and not in SPOWER 
(which computes P/C) the following approach will be followed to determine 


if the sun is visible before the force is computed. 


First, the angle between the 
relative sun and the radius 
vector is 


B= cost] -¥ . (BF 
r Ys- Ty 
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Then, the angle defining the limits of the cylindrical shadow can be 
\ defined as 


3 6: sin™t (2s) 
‘ r 


and the requirement for visibility is 


43> O 


There are, of course, several assumptions which have been made which, 
though not completely apparent, should be noted since they introduce 
slight errors 


1) The sun is assumed to be a point 

2) Ibe earth is assumed spherical for the purposes of 
defining the shadow 

3) The shadow is asswzed cylindrical 


These assumptions are believed to be reasonable for most applications 
and should introduce negligible errors due to the small effect of Fs 
on most trajectories, 
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Computational Logic: 





Position of sun relative 
to satellite 


=r 





=~ —_ 
R. = r 


Check for shadow 


coss= r-(-R,) /r Rg 
sin =/1-c0s28 
& =tan-L sin | 
COS 





sing =R,/r 


cos 6 =\VJl-sin 6 


8 =tan7} (ss 8‘ 


cos 8 
@ 743 <> 
XS = 0 K 
F. 
RETURN 
END 
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Then, the angle defining the limits of the cylindrical shadow can be 
defined as 
©: sin"+ (2s 
ry 


and the requirement for visibility is 


7 >89 


There are, of course, several assumptions which have been made which, 
though not completely apparent, should be noted since they introduce 
slight errors 


1) The sun is assumed to be a point 

2) The earth is assumed sphericel for the purposes of 
defining the shadow 

3) fhe shadow is assumed cylindrical 


These assumptions are believed to be reasonable for most applications 
and should introduce negligible errors due to the small effect of F, 
on most trajectories. 
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Computational Logic: 


Position of sun relative 
to satellite 






Check for shadow 


cos@= ¥-(-R) /r R, 
sin = JV 1-cos28 
& =tan-1 sin < | 
cos 4 


sino =R,/r 


cos9 =VJl-sin © 


8 =tant ‘ee 2) 
cos 8 





278 os 





RETURN 
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2301-24-11 Function SPOWER 


Purpose: . WER computes the ratios of the solar power to the 
vzed of light in mks units. 
vo allow for .otion about bodies other than the earth 


Deck Name: SPOWER 


Calling Sequence: SPOWER (RS) 


Input/Output: 
FORTRAN 
I/O Name 
T RS 
0 SPOWER 


Subroutine Required: 
Functions Required: 


Approximate Deck 
Le~gth; 


AE AS TS 


This routine was intended 


Math Common/ 
Name | Dimension Argument Definition 
rg-T 3 Arg Position of sun 
relative to the 
satellite 
P/c 1 - Ratio of solar power 
to the speed of light 
None 
None 
50 (octal) 
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2436 keee5 Subroutine PERT 


Purpose: PFRT is 4 special rurpos® routine which compytes the 
gravitational accelerations (perturbative) af the sur 
and moon on a geocentric satellite, The assurption is 
made in the develonment that the seorentric radial tis- 
tance is small compared to *he distance +n the disturbing 
mass: thus a more exact formulation is required for the 
general case, 


Deck Name: PFRT 


Calling Sequence: Call PERT (ACCFI,) 





Tnout/Ousput: 
FORTRAN Math Comnon/ 
T/O Name Name Dimension Argument Definition 





R T 3 WRK (44) radivs vector in the 
: frame of 1950.0 
RS Ts 3 WRK (98) radius vector of the 
4 sun (moon) in the 
‘ RM r 3 WRK (101) frame of 1950,0 rela- 
® tive to the Farth 
TM An ] CON (&) gravitational constant 
US Mas ve CON (9) for the moon and sun 
cs ad ‘ i“ 
ACCEL at 3 Arg perturbative accelera- 
tion vector due to 
gravitational accel- 
eration of the sun and 
moon 
Subroutines Required: None 
Functions Required: AMAG (vector magnitude), FQ (Fneke series Minction) 
Approximate Deck . 
Length: 197° (octal) 
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Formulation: 


The acceleration experienced by a body (traveling about the Earth, for 
example) resulting from masses other than the central body (i.e., the sun, and 
moon) is the summation of terms of the following form, 


a «|2-8,-8 


Fr Ri}? i 


This equation, while correct, exhibits a very undesirable nature for 7<<¢& 
since for these cases the term in the brackets is the difference of two nearly 
equal quantities, Therefore, since this limiting case is the case of primary 
interest, an alternate form is necessary if these terms in the acceleration are 
to be included. One zuch modification is based on the similarity of this per- 
turbatien equation and that for the Encke acceleration. 


consider: 


5 ee 7 = 7 = 
r 4[-4 +B) Fal R-r 


elP)| 


“< [f a Ge 9)? | 


gl 


‘RS 


where ag before, since f£g has one more degree of freedom than the ratio Rlo 
g can be defined as follows 


I- fg =(1 +2) * = Ros 


or k 
sche. /- 
fea CLS |) 
= 3p- 3:5 9* eS. 3 FP -aiiy 7 ¥ ee 
ZG * Sg’ ten g' + 
also 


a= /t83 
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2 ~ gps (R***) 
but 
7 = R-F 
= (X-x)X--- 
. therefore 


| ge wah (Wi P+ B*(X-xY (Fy) -@-z)\) 


“ape (2( ar y #22) “(x8 ry? 24) 


‘ Computational Logic: 













| eet ae 
compute poSitlion of sun and moon 
relative to the vehicle 


compute Encke series for solar and 
lunar perturbative accelerations 


construct resultant perturbative 
acceleration 
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246301-2-5.1 Function FQ 


Purpose: to evaluate the Encke series utilized in constructing solar 
and lunar gravitational accelerations of close geocentric 
satellites, 

Deck Name: FQ 


Calling Sequence: FQ (RN, R) 


Input/Output: 
; FORTRAN Math Common/ | 
I/O Name Name _| Dimension | Argument Definition 
a RN R 3 Arg radius vector for the 
sun (or moon) in the 
frama of 1950,6 
I R F 3 Arg radius vecter for the 
vehicle in the frame of 
1950,0 
0 FQ Le 1l ~ Numerical value of 


Encke series 


Subroutines Required: None 
Functions Required: None 


Approximate Deck 
Length: 120 (octal) 
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2.3.1.2.6 Subroutine EPHEM 


Purpose: EPHEM is designed to compute approximate position 
vectors for the sur and moon from two-hod:r eqnatiors 
of motion and tabrtated position and velocity Inform 
ation obtained from the Data Input Group. 


Deck Name: EPHEM 


Calling Sequence: Call EPHEM 


Tnput/Output: 
FORTRAN Math | Common/ 
T/o | Name | Name | Dimension Argument Definition 
T TW t 1 WRK (50) The whole and functiona? 
part of a mean solar 
TF 3 WRK (51) day defining the desired 
epech for F, and ff, 
0 RS t 3 WRK (98) Solar and lurar posi- 
*ior vectors in the 
RM Tr 3 WRX (101) | coordinate frare of 
1950.0 
T GME, yo 7] con (4) gravitational corstants 
GMM Pin zd con (8) for the Marth, the 
GMS Ps st CON (9) moon ang the sun (in 
Km3/sec“) 
if FPHAM Ty Fy 339 FPHOM Tabulated position 


and velocity vectors 
for the sun and monn 
for the time periad 
1965 — 1975 






FPHAM (L) Subarray contaiving 


the Junar data 





EPHAM (274) | Snbarray cortaining 


the solar data 
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Subroutines Required: CONIC (conic section routine) 
Functions Required: None 


Approximate Deck 
Length: 456 (octal) 


Formulation: 


The equations of two-body motion (i.e., the Earth-sun and the Earth-moon 
systems) relative to the center of mass are 


pies a 2 
i Ging (72. 
mrm, 


“By 


BN ” a ~— 
ly =-Gm (Aa ig 
"(MN 447, 3 

n 


which are immediately recognizable as the equations of conic motion, This fact 
will be utilized in the following paragraphs to approximate the true ephemeris 
of the moon and sun relative to the Earth. 

Consider a table of values for *% and Mat regular intervals of 
approximately 1 year from which the conic describing the motion of the sun 
(relative to barycenter) can be computed. The problem is simply to locate the 
closest tabulated point (no more than about 5 months), define the time over 
which the body has moved, and compute + and V utilizing the conic section 
routine developed for the reference trajectory, Since this conic is altered 
but a very little by planetary perturbations, the result is of moderate Q 
precision. 


Now consider the problem of the moon's motion when the sun produces a 
significant perturbation, First, it is obvious that the time interval between 
values of T, and ¥% must be reduced so that the effect of the perturbations 
can be kept small (for the present case At = 3 months). Second, it is apparent 
that the results obtained with the technique would be more accurate if they 
could include some of the effects predicted if the next-to-closest point in the 
table were utilized for extrapolation purposes, This correction was included 
by weighting two separate computations of the position vector (in a manner 
which resulted in even weighting for a date midway from either point in the 
table and a weight of 1 for the term designed to include the effects of the 
closest point if the time was equal to a tabular time) and averasineg the 
results, 


The final step in both cases was to translate the origin from the bary- 
center to the geocenter, 


Results prepared by this routine have been checked against a true 
ephemeris to compute the errors inherent in the process. This analysis showed 
that the sclar position data ware completely satisfactory (from the standpoint 
of accuracy) for the computation of both the solar radiation force and the 
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solar gravitational acceleration in that the noise in the data would normally 
produce errors in the total perturbating acceleration (from MOTION) which were 
beyond the 8-digit limitations of the machine, The same is not quite true (for 
near-Earth satellites) of the lunar ephemeris even though the correction cvcle 
is employed, However, it is noted, in defense of its use, that the second 
intepral of the total resultant differential acceleration ( 4r ) is in itself 
a correction to the conic position vector ( % ). Therefore, the noise which 
is not lost when summing the perturbative terms in MOTION is generally lost 
when adding ar to % . 
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locate next’ 
closest point in 
table and store 

lunar position vector 


edith at OE sy ee et ee 


locate closest 

tabulated value of 

position and velocity 

for the sun from data 
stored in common. 

Then using the 2 boa, 
equations compute the 
position of the sun 
relative to the barycenter 


repeat the process for 
the moon 


<> 
2 values of position 
are avaliable 


compute weighted average 


of two lunar positions 


transform all data from 
the barycenter to the 
geocenter 
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2.3.1.3 THE JNTEGRATION GROUP 


The differential acceleration vector (relative to the central force 
acceleration describing the reference ellipse) is integrated to define 
poth the differential position and velocity vectors with a Gauss-Jackson 
backwards difference formulation. This solution is accomplished by 
mechanizing several separate routines. 


INGRAT The driver routine 

HSTZE Dynamically varies the step size for the 
integration to maximize accuracy and efficiency 

START A fourth order Runge-Kutta integration routine 
utilized to establish the required difference 
table. 

DIFTAB Evaluates the leading diagonal in the sum- 
difference tables 

INTEG Gauss-Jackson backwards difference uncorrected 
integration 


The various portions of this group are in themselves general purpose 
routines. However, the group as a whcle has been designed as a special 
purpose routine for the differential corrections program. This design 
philosophy can be understood when the requirements for the routines operation, 
the nature of the acceleration being integrated, and the requirements for 
communication between the calling program and INGRAT are analyzed. This 
approach also affords the maximum efficiency from this integration concept 
by eliminating checks for alternate stops, etc. , 
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Subroutine INGRAT 


Purpose: INGRAT serves as the driver for the integration of the 
equations of motion employed in this program and is 
intended to function in the mode of a Gauss—Jackson 
integration. Complete logic is included for starting 
the integration and/or varying the step size and re- 
starting the process, 


Deck Name: INGRT 
Calling Sequence: Call INGRAT (ISTART) 


Input/Output: 


FORTRAN Math Common/ 
I/0 Name Name | Dimension Argument Definition 


t/o ISTART ~ 1 ARG A fixed point index 
utilized to define 
the mode of operation 
for INGRAT. ISTART = 1] 
| for Runge-Kutta starter; 
= 2 for Garss—Jackson 


continuation 
if CONV1 sec a) CON (11) Conversion fren mean 
day solar days to seconds 
a RCO Te 3 WRK (28) Position and velocity 
vectors for the conic 
vco ¥ 3 WRK (31) reference trajectory 


at the time of the 
most recent rectifi- 
cation (Kn,Km/sec) 

in the frame of 1950.0 











t/o WRK (LL) The position and velocity 
vectors in the mean 
WRK (47) equator of date frame 
of 1950.0 (Km, Km/sec) 
t/o WRK (50) The whole and fractional 





number of days elapsed 
since the reference date 
1950.0 (JD 2433282.423) 
for the present epoch 





WRK (51) 
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Common/ 

Argement | Definition 
WRK (58) The acceleration, 
WRK (55) velocity, and position 
WRK (52) vectors defining the 
true motion relative 

to the reference ellipse 


in the frame of 1950.0 
(Km/sec®, Km/sec, Km) 





The whole and fractional 
number of days since 
1950.0 for the data 




















The step size in sec~ 
onds for the numerical 
integration 





The position vectors 
for the sun and moon 
in the form of 1950.09 
(Km) 





The time interval in 
seconds since the refer- 
ence trajectory was last 
rectified 


Subroutines Required: CONIC (conic reference trajectory) 
START (Runge-Kutta starter) 
DIFTAB (difference table) 
INTEG (Gauss-Jackson integration) 
HSIZE (step size control) 


Functions Required: AMAG (vector magnitude) 
Approximate Deck 

Length: 501 (octal) 

Discussion: 


INGRAT is a special-purpose driver routine for a step-by-step numerical 
integration of the equations of motion as developed for this program. This 
integration ig accomplished by mechanizing a step size check routine and either 
a fourth-order Runge-Kutta starter or a Gauss—Jackson (aouble-sum) backwards 
difference (inrough sixth difference) contimation routine, 
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INGRAT accomplishes its objective by continuously monitoring the time 
interval f:0m the present epoch to the time at which data will be available 
and defining the information required during the checks performed within HSIZE 
to assure that an optimum step size can be utilized to arrive at the target 
epoch exactly in the fewest number of Gauss-Jackson steps possible. In addi- 
tion, INGRAT serves as the means by which data generated in START (the fourth 
order Runge~.utta starter) are stored and differenced to set up the solution 
by the Gauss-Jackson routine, 


Several steps are taken to assure that the solution proceeds in an 
accurate and efficient manner, However, no corrector cycle has been employed 
either in the starting process or in the continuation process. This step is 
felt justifiable since: 


l. The equations being integrated represent the displacement from and veloc- 
ity relative to the reference ellipse (Thus, even moderate errors will not 
generally be realized when the true position and velocity vectors are com- 
puted by adding the results of the integration to the conic position ar 
yelocity.). 


2, The Gauss-Jackson formulation has demonstrated superior performance capa~ 
bilities when compared to the Runge-Kutta formulation in the integration 
of X = -KX for similar step sizes without employing the corrector. (For 
this reason, the Runge-Kutta method [which is employed to start the Gauss- 
Jackson integration] will utilize a step which is one-half of that to be 
employed in the Gauss-Jackson method, ) 


3. The step size itself will be controlled in such a manner that the predictor 
formula alone will provide the required accuracy by monitoring the contri- 
bution of the last term carried in the integration series to the integral 
(this procedure will cause the step size for this routine to be less than 
that for a predictor-corrector formulation; however, sample calculations 
have shown that the slight improvement in accuracy [about one digit] and 
step size [a factor of about 2} ) does not justify this procedure for this 
application, 
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DEL6* MAxX(AS) 
DEL S* MAX(A’) 
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CALL HSIZE, 
ISTART= 2 
r, &/ ; ART 
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LIMIT > | 

IENOF 2 






HRK = ‘Y2 
HH ® WRK 


CALL Conic 





rxy=/r-rcHack] 

Ts TtHRK Alias 
; PVvéC 2 /COMICer 
CHECK F! ICHECK */ en lr sroresr WE = WCONIIC #°§ 








ICHECK =2 






LIMIT #7 


STORE(Limur, 1) = PI) 
LUMIT = LMT +4 


ce Limreg 


CALRK DIFTAB 


PCHECK 7Pr 
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263e1.3e1 Subroutine HSIZE 


Purpose! HSIZE is a step-size, control routine designed to operate 
in conjunction with the Gauss-Jackson integration pack- 
age. In its simplest form, HSIZE determines whether the 
stepsize (h) should he halved or whether it can be doubled 


without loss of accuracy or chance of producing a loop. 


Deck Name: SIZE 
a Calling Sequence: HSIZE (TTOGO, DEL6, ACCEL, RXX, ISTART, L, IH, H) 
Input/Output: 


FORTRAN Math Common / 
I/0 Name Name Dimension | Argument 


a TTOGO t x ARG 
go 


DEL6 nf ARG 
ACCEL z ARG 
RXX BS ARG 





ISTART ARG 
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Definition 


the time interval (sec) 
between the present 
epoch and the epoch of 
the next piece of 
observed data. 


the magnitude of the 
largest component of 
the vector composed of 
the sixth differences 
in tne acceleration 


(AT = Km/sec2). 


the magnitude of the 
acceleration vector 
(Km/sec2) 


the magnitude of the 
change in the second 
integral compared to 
the previous step. (Km) 


an index used to define 
if INGRAT is functioning 
in the start mode 
(ISTART=1 for Runge- 
Kutta) or the continue 
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Common / 
Argument 

































Dimension Definition 


FORTRAN Math 
Name Name 





mode (ISTART=2 for Gauss~- 
Jackson) 


an index used to check 

| the number of consecutive 
times through HSIZE that 
the step size has been 
adjudged to be too small 
(when L=4 the step is 
doubied ) 


an index used to deter- 
mene if a step is to be 
| guessed (to be refined 
later) or if the present 
step is to be checked, 
TH=1 for guess; =2 for 
present step 


the step size to be 
employad in the numerical 
integration of the equa- 
tions of motion employing 
the Gauss-Jackson back- | 


wards difference formula- 
tion, (sec) 


Subroutines Required: None 
Functions Required: ALOGLO (log to the base 10) 


Approximate Deck 
Length: 377 (octal) 


Diecussion: 


The basic rationale for deciding whether the stepsize for the integration 
(Gauss) is adequate is based on a set of more or less arbitrary limits for the 
absolute value of the largest sixth difference (DEL6). Thase limits, in turn, 
are based on the order of the magnitude of the acceleration and were selected 
in such a manner that the integration series would be convergent to the desired 
degree with terms through the sixth difference included, 


One slight problem arises in the application of this rationale in that if 
the acceleration itsalf ever passes through zero, the tolerance goes to zero; 
and the stepsize is immediately adjudged to be too large. To surmount this 





problem, a second test is employed in those cases where the quantity DEL6 is 
larger than that allowed by the upper limit. This second test requires that 
the contr shubton of the sixth difference to the second integral be less than 
one part in 104 of the change in the second integral for the last step, i.e. 


6h 
Ar ®& Hae / 15 


< a r on 
: | Ge Ata | 
or say 
3 2.6 > “a 
oa’ | ats - at, | 


(The seccnd integral is the displacement from the reference conic which is 
growing in a reasonably slow fashion.) The step is always halved if the upper 
tolerance level is violated. 


Three other modifications to this basic rationale were also coded for the 
sake of operation, efficiency and accuracy, Each of the threo pertains to the 
case where a signal would normally be generated to double the step size (DELS 
less than lower limit), The first modification is employed to assure that the 
step size is too small at several consecutive points before a signal is 
generated to double the step, This practice assures that fewer cases in which 
the tolerance range is temporarily violated on the low end will result in 
restart. The second modification is employed to assure that after a whole 
number of the new integration steps, the time-to-go can be reduced to zero 
(this performance is essential in the operation of the program). This second 
test is passed automatically if the step is halved; however, if the step is 
douvled, time-to-go divided by the old step must be even, The third modifica- 
tion is employed to assure that the integration can be effectively restarted 
prior to the time that the time-to-go will reach zero. 
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2.3.1.3.2 Subroutine START 


Purpose: 
Deck Name: 
Calling Sequenc 
Input/Output: 
FORTRAN 
| 1/0 | Name 
I RCO 
vco 
) RAD 
VEL 
1/0 RDD 
RD 
R 
{ 
1/0 TIME 


e: 


START is a fourth order Runge-Kutta integration routine 
utilized to establish the accelerations at the seven 
points along the trajectory and define the difference 
table required by the Gauss-Jackson integration. 


START 

Call START 

“Math ~ 

Name | Dimension 
Teo 3 
Vo ° 3 
50 3 
¥50 3 
ar 3 
Pe 3 
ar 3 
t-t, 1 


198 


. 


Common/ 
Argument 


WRK (28) 


WRK (32) 


WRK (4d) 
WRK (47) 


WRK (58) 
WRK (55) 
| WRK (52) 


Definition 


The position and vel- 
ocity vectors in the 
frame of 1950.0 des~ 
cribing the conic 
reference trajectory 
at the time of the last 
rectification (Km, 
Xm/sec) 


the position and velo- 
city vectors in the 
computational coordinate 
frame of 1950.0 

(Km, Km/sec) at the 
epoch tth 


the acceleration 
(Km/sec*) the velocity 
(Km/sec), and the posi- 
tion (Km) vectors 
relative to the refer- 
ence ellipse in the 
coordinate frame of 


the time interval in 
seconds since the last 


rectification of the 
reference ellipse 
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CONI 
MOTI 


Subroutines Required: 


“T Comnon/ 
Argument 


WRK (117) 


Definition 


the instantaneous 
reference positior 
and velocity vectors 
(1950.0) (Km,Km/sec) 


WRK (120) 


WRK (123) the step size employed 

in the integration 

(sec). This quantity 

is one~half that employed 
in the Gaussian inte- 


gration , 


C (reference conic) 


ON (total acceleration vector relative to ellip- 


tical reference) 
Functions Required: None 
Approximate Deck 
Length: 500 (octal) 
Discussion: 


START functions in the 
That is, the second-order d 


standard single integration Runge-Kutta mode. 
ifferential equation is reduced to two first- 


order differential equations with the substitution 


bh 
<abay: 


and each part which may be 
x= 

as integrated as follows: 
x (t 


where: 


i es 


f2 


i 


£3 
fi, 


- 


represented as 


PCG) 


t+h)= x (t) + b (£1 #2f, + 2f, +f, ) 


£4%,%) 
f (tth, xth fz) 
2 2 


f (tth, xth fo) 
ae 
f (tth, xii f3) 
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No corrector cycles or smoothing operations are superimposed on this pro- 
cedure, Thus, Since the errors tend to grow more rapidiy in this method of 
integration than is the case with the Gauss-Jackson integration, it was decided 
that the step size utilized for each step of the starter would be no more than 
one-half of that utilized for the continuation integration. This step (made 
in INGRAT) assures a more consistent level of accuracy in the two routines. It 
is noted, however, that the subroutine MOTION must now be called eighs times in 
this mode of operatioi, compared t+ one time in the Gauss mode, 


A complete discussion of the development of the Runge-Kutta family of 
integration formulae can be found in many texts treating the subject of 
numerical integration, One such text is Numerical Analysis by K, Z, Kunz and 
published by McGraw-Hill in 1957, ee ge 
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6a 
T= TIME +H/z 
OTT 


CALL CONIC 


= F CONIC 
= VCOAIIC 
=H® RD 
2H*® ROD 
XK *RtFI/2 
Xo=R0+G//2 
YFAZRCHY 

VE =VC#XO 

T= TIME +H/2 
O1?T 






EVALUATE CONIC 
Po AND V AT r+H/2 

















F2 7H XD 
G2 =H*XDD 
X*R4F2/2 
XDO=RDtGA/2 
RA ROX 
VE = VC+tXD 


am) 











EVALVATE CONIC 
— eed 
ro ANDO VY AT TH 









F3 =H# XD 
G3 = Hx* Xb0 








x2 ARTF3 
XD3* RO#G3 
RA= RC #X 





VE= yVOrxX 
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F4Y=H#XD 
GY 2H* XDD 
RERtFY/6 tF 2/3 tFY/3+FY/E 
RD=RDtG//6 + G2/3+G3/3 +G4/5 
RAD= RCA/ +R 
VEL= VCA/#RD 


RETURN 
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263e143.3 Subroutine DIFTAB 


Purpose: 


Deck Name: 
Calling Sequence: 


Input/Output: 


RSTORE 
VSTORE 





DIFTAB differences seven values of the acceleration 
vector and constructs the diagonal of backward differences 
and leading sums required in the Gauss-Jackson integration 
routine (INTEG). 


DIFTAB 


CALL DIFTAB (R°TORE, VSTORE, STORE, H, COEF) 


Common/ 
Dimension | Argument Definition 


the position and velocity vec- 
tors relative to the reference 
ellipse at the center point in 
the table. (Km, Km/sec) 


the array of acceleration 
vectors at seven times (spaced 
at intervals of h) which are to 
be differenced and summed 
(Km/see?) 


the step in time which correlates 
the various entries in the STORE 
array (sec) 


the array of numbers which is 
required by INTEG for the inte- 
eration of the equations of 
Motion. The first subscript 
denotes the second sum, the 
first sum, the first difference 


the sixth difference. The 


second subscript denotes the 
xX, y and z components. 


eeery waane Sixuo 


Subroutines Required: None 


Functions Required: 


Approximate Deck 
Length: 


None 


470 (octal) 


208 
S™ 65-1203-3 


ea se Sal erttieatinet: Aeris shal Alli» aeeen naan ccopeetiinie pci nan a a Fa a a pA ne 


Discussion: 


The output of the start cycle is an array of acceleration vectors at 
seven times (spaced h seccnds apart) and position and velocity vectors for 
the center point. These numbers will be utilized to construct the leading 
diagonal requir>4 in the Gauss-Jackson integration process as follows 


1) The six differences of the components of the acceleration vector 
will be constructed. 


2) The second sum corresponding to the position for the central value 
of the acceleration is calculated from: 





oe 2 
OMG Re Be. Bike 
ke a a 
: ° h* 12 240 = 60480 


3) The average value of the first sum ccrresponding to the velocity for 
the central value of the acceleration is celculated from: 
ieee aon + Ay) (Ary, +3 1A» +A¥%) 
°o fh ZY 1440 
F 190(A%, + A$) 
120960 


and the corresponding first sum on the proper diagonal in the difference 
table is 


4) The first and second sums are stepped down to the next to last 
» diagonal and the terms along this diagcnal are loaded into an array 
set aside for input into INTRG. (The next to lest Qiaronal was selected 
since INTEG undates the diagonal before performing the integration). 


In order to preserve accuracy in these operations the differencing will 


be performed in double precision. This procedure has been shown to conserve 
approximately one significant figure when opereting with the IBM 7094. 
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Computational Logic: 
















CONSTRUCT DIFFERENCE TABLE 
FOR EACH COMPONENT OF 





USING THE POSITION AND 
VELOCITY VECTORS AT THE 
CENTER POINT IN THE TABIE, 
EVALUATE 3’, ©“ AND STEP 
THESE SUMS DOWN TO THE 
DIAGONAL OF THE NEXT TO LAST 
ACCELERATION. THE GAUSS 
CONTINUATION WILL STEP ALL 
NUMBERS TO THE NEXT DIAGONAL 
BEFORE STEPPING THE FUNCTION 


EQUATE THE COEFFICIENT ARRAY 
(DIMENSIONED 8 BY 3)TO THE 
DIAGONAL CONTAINING THE NEXT 
TO LAST VALUE OF X, Y, % (THIS 
VECTOR CORRESPONDS TO R, V AT 
THE LAST STEP TAKEN BY THE 
STARTER ROUTINE ) | 
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2630146344 Subroutine INTEG 


Purpose: INTEG continues 2 stepwise intogration of 
the equations of motion in a Gauss-Jackson 
mode once started by an independenw process. 
The routine employes six diffences and is 
formulated around the backwards difference 
concept (no corrections) 


Deck Name INTG 
Calling Sequence: JALL INTEG (SAVE) 
Input/Output 
















WRK(71) |The array containing 
the diagonal to be 
utilized in the next 
integration step em- 
ploying Gauss! Equation 
(2 sums, 6 differences, 
3 components of 
acceleration) 








WRK(58) |The acceleration vector 
defining the motion 
relative to the refer- 


ence ellipse (Km/sec@) 

















i/o SAVE Ar, Arg |The acceleration vector 
at the last By eerste 
| step (Km/sec“) 
0 RD Ar WRK(52) |The incremented 
R Ar WRK(58) |position and velocity 


vectors for the motion 
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FGRTRAN Math Common / 
Name Name Dimension Argument 









Definition 


The stepsize for the 
numerical integration of 
the aq. of motion (sec) 










————— - 


STEP 


Subroutines Required: 


_ 


None 





Functions Required: None 


Approximate Deck Length: 330 (octal) 


Discussion 


The continuation of the intregration will be performed by employing 
This process is accomplished 


a Gauss-—Jackson or double sum formulation. 
by mechanizing the following backwards difference fcrmulae 





e é ee 2 
. 8 , tes 3A. 
ary H(Z,y, "ee easy - gus 
+251 AS +95 AY 4+ 19087 Ars 
720 '-% 288 (2 60490 *° %% 
+5257 Ae 
17280. ‘3 
eT Oe eee eee 
we oT TB TZ! “agg 
3 4 ; s 
*la2 Qk. +/726. A’ +45 Ae 
ayo °° Begg "ie R492 t-% 
+ 1s. 2° ) 
2¥0 «+3 
i a ax 
Cth t-% 
2 2 ? 
wal es 
cee = L ¢ Le 
i ae oo 
Ais = x; 7 Key 
2 s U t 
A;., Avare 43. 
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This mechanization will be aided by the establishment of a two dimensional 
array which will be referred to a8 the coéfficient array (COEF) and which 
will contain the information required by the formula. 


M 


Zz 
COEF (1,1) COEF (I, 2) * diy Boats 


COEF(2,!) = Ly 


COEF(3,)= AY 


COEF (8,/) * Ay 


(This array is thus seen to contain the terms along the second leading diagonal) 
Upon entry to INTEG, the diagonal which is stored in COEF will correspond 

to the last integration step; thus, the first task to be performed will be 

the updating of COEF according to the definition of its members. This 

done, the trajectory can be stepped. 


Two specific comments are required regarding this process. First, 
the process is obviously not exact and thus for extreme accuracy a 
corrector cycle based on central differcences should be employed. This 
type of operation has been avoided by attention to detail in the formulation 
of the trajectory portion of the prorram and by careful selection of the 
stepsize. However, this revision is felt tc be advisable if any numerical 
problems develop during the application of the differential corrections 
proprram to the determination of satellite orbits. 


The other comment is an observation which pertaines to the fact 
that two integrations which are performed are both based on the acceleration 
rather than on the acceleration and velocit;. This fact results in the 
evaluation of the second integral to an accuracy superior to that obtained 
with the other approach or for that matter, to that obtained for the first 
interral. This behavior is vers desirable for this program since it assures 
that the major contributives to the acceleration at future points along the 
trajectory (primarily functions of position) will be well known. Care must 
be taken, however, in applying this logic to trajeccories for which the 
non-linear effects of drag, etc. are involved since for these cases the 
accumulative errors in the first integral (X) could drastically effect the 
predicted path. 
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Computaticnal Logic: 

















STEP THE TERMS ALONG 
THE OPERATING DIAGONAL 

IN THE DIFFERENCE 

TABIE SC THAT THE PRESENT 
VAIUE OF THE ACCELERATION 
IS CONTAINED. THEN SAVE 
THE ACCEIERATION FOR FUTURE 
STEPS. 


EXTRAPOLATS THE MOTION 
| USING GAUSS! FORMULAE 
FOR THE FIRST AND SECOND 

INTRGRAIS. SIX DIFFERENCES 
CF y ARE CARRIED 
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2.3.1.4 THE TRACKING GROUP 


This group of routines is desigried to define the position and velocity 
of the satellite relative to each of a set of prescribed tracking stations. 
This task is performed by defining the position vectors for the group of 
tracking stations as a function of the time, and associating with each, 

a set of unit vectors describing a topodetic coordinate system. This 
information is then utilized to determine the range, range-rate azimuth 
and elevation of the satellite relative to each of the stations. 


The operations are performed with the aid of the following routines: 


TRAK Driver routine for computing the relative 
position and velocity data for each station 


GHA Defines the position of the Greenwich meridian 
as a function of universal time 

UNIT Constructs the position vector for the 
tracking station and the topodetic coordinate 
system 

EQINOX Computes the correction to the computational 


coordinate frame of 1950.0 (used for trajectory 
definition) to adjust for nutation and 
precession 


Following the construction of the relative position information for 
a given station, a check is made to determine whether observation data 
are available for the epoch under analysis for that station. If not, 
the checking of the various stations is continued. However, if data 
are available, transfer is made from this group into the Filter group 
and a correction to the state vector is computed before continuing the 
station check process. This link between the two groups is very 
important and deserves considerable attention. 
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Purrose: 


TNerlk Mame: 


“Salline Sequence: 


Trent /Outoyt: 


ra I 
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rr 


FORTRAN 


Nama 





TSTAR” 


CONT. 
nana 


MONT 


vourcy 


HONDO 


sirarn 


a ee 


Sukrontine TRAX 


eee 


™IAY is desizned to checlh each of the trachinz stations 
brinzt errleved at each noint alorz the tratector: for *he 
rurpose of idextifving those stations at which the vehicle 
is visible and to check the data tane to see if data is 
avaitable at that tire before transferrine to the data 


filter, 
PRACY 


Cali TRAY (TSTART) 


Math 
“are 
- a Art 
es 1 cov 
= 7 nor 
= 7 COM 
= ] SAT 
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Cormon/ 


Armmert 


(12) 
(12) 


(1h) 


(16) 





Index utilized to 
restart the intezration 
when the next data noint, 
has been read irto 

r ano Mm 


Conversion factors 
(davs to seconds and 
degrees to radians) 


Ojtrut tare number 


Trdex utilized to check 
stations for visibility 
at al? times or just at 
ths times for which there 
is data 


Tndex to limit number 
of stations cneced 


Array containinz data 
for the stations to ba 
checked (latitude, 
loneitide, altitude 
and nane) 





STN 65-1203-1 


ae 


1/0 


T/9 


FORTRAY Math Gorron/ 
Name Nare | Mimension Arcunert, Hsefiri tion 


HORCOR AFI 19 cay aa Gra Herizor corrections 
“or “ach of the ur to 
1€ stations ers lorcad 


Number = 2 STA (241) “otal mumbcr o% stations 
bein? erriored 


ROTATE NP 3 x 3 wex (1) | Matrix for trans *ormirs 

ROTINV PINT 3°x.3 WRY (10) tre frame of 1950,C to 
that of data and its 
inverse, Th = PPP en 

EN N 3x3 WRK (19) Mutation array relating 
the mean equator of date 
to the tre equator of 
date, Th = YY, 

TTRANW ter 1 WRY (12) Time of jast data roint 

TTRANE ] WRK (133 

RVFC, Ts50) 3 WRK (hb) The cartesian resition 

VES V5 3 WRK (47) and velocity vectors of 
the satellite at the time 
TW“ TF in the reference 
frame of 195C,° 

TW, TF % digs. WRK (50) The two words d°firinz 
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Tracking station lati- 
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(arc tangent) 
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Description & Formulation: 


TPAK is designed to check each of the stations being employed in the 
tracking network (one to ten) to determine the relative position of the 
vehicle at that epoch. As a portion of the process the range, range-rate, 
azimuth and elevation of the vehicle in topocentric coordinates are 
computed for the purpose of comparison with the actual tracking data in 
a differential correction of the trajectory itself. The procedure 
involved is presented below: 


Upon defining the position and velocity vectors in the frame of date 
and computing the Greenwich Hour Angle associated with the time (GHA), the 
locatior of each of the stations is computed and a set of unit vectors 


defining a topocentric (radar Az-EZ) coordinate system is constructed (UNIT). 


At this point the relative position and velocity are computed 


Arto ep el 4 f= goly-tp 


and the elevation and azimuth are defined as follows: 


sin E2 = P ‘U -30< E£<90 
Pp 


cos EQ = + 1-sinZEg 





aA 

cos a = f-N = cos A, cos E2 
P 
> A 

sin b sin A, cos ER 


it 
~> 
m 

tt 
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tan A, = Pp: 


Wy 
ee 


When the data point has been reached (i.e., the time at which data 
from the tracking stations is available) and associated with the proper 
Station, the FILTER group is called and the data point processed, However, 
before the solution can proceed and before the trajectory can be stepped, 
the radius and velocity vectors in the frame of 1950.0 must be defined. 
This step is necessary due to the fact that subroutine UPSTAT may have 
adjusted the trajectory discretely in order to correct the tendency for 
the observed and computed trajectories to diverge, 


TRAK is designed in such a manner that the routine can be by-passed 
in two distinctly different modes, First, if specified in the input data, 
the stations can be checked at every point along the trajectory (i,e., at 
every integration step) so as to provide a record of the track or only at 
those times when data is known to be available. This option drastically 
improves the general efficiency of the computational logic, Secondly, if 
desired, a minor efficiency can be effected in that only that staticn 
actually recording data at the time of the observation need be checked. 
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, I = ITRACK 
CALL STMAT 


Trogo = TpaTa-T 
TSTART = 1. 


1 = (Taonga /H) + 1 
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2,3.1.4.1 Subroutine UNIT 





Puyroge: WNTT corrntes the rnsition vector defirine the 
Taratian of the +rackinz stations and evalnates the 
comments of the (im, Fast, North) unit vectors at 


the station, 





Neck Mamo : UNTT 
Callins Seqierce: Call INT? (CHA) 
Tart /Ontont 

BOPTRAN 
T /0 Mame Nimersion 
T GHA GHA 1 
7 RF 2. 7 
T RPOT Rp 1 

2 

(a) RT rp 3 
tT STAT Ty ) 
T STON X . 
F SAT™ H 7 





2038) = 


~ 


Wry 





(195) 


(104% 


(197) 


Nefinition 


Hour Ancte of Creenwich 
relative to tyme yverng] 
~qiinox (rad) 


Sqvatorial radivs of 
the earth (Yr) 


Palar radius of the 
earth (Km) 


Positjor yertor Por the 
trackine station at 
this time (%m) 


Station Jatitrde 
(Cendetic) ir radians 


Station loatitude in 
radians 


Station altitnde (%m} 


Yn, Bast, North unit. 
vectors exnressed in 
nartesian coordinates 
in the tme equator 
of date frame 
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Subroutines Required: None 


Functions Required: SIN (sine) 
COS (cosine) 
SQRT (square root) 


Approximate Deck Length: 170 (octal ) 


Formulation. 


UNIT computes the components of a set of topocentric unit vectors 
for each of the tracking stations being checked as a function of the 
time and the radius vector of the tracking station itself. The unit 
vectors will be constructed first utilizing the information and defi- 
nitions presented on the following sketches. 
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ae NI ree nee SE Re EE renner rome gemene reaper mit me 
‘daa tt op ees = — 








thus 


U cos L O sin L cos R.A. sin R.A. O of? 
E] =& Q i (a) ~sin R.A. cos R.A. O Y 
N -sin L O cos L Ye) fe) ' nN‘ 


cos L cos R.A. cos L sin R.A, sin L Yi 
-sin R.A. cos R.A. 0 Y 
-Sin L cos R.A. -sin L sin R.A. cos L N 


ut 


The vector defining the position of the tracking station is, however, 
more difficult to construct since it is not found by simple rotation. 
This construction is, however, simplified by noting that 

os > A 
YS 7 or Wy 
: A 
= r(cos 8 Ro + sin & N ) + nu 


and that the quantity f can be evaluated by considering the equation of 
the ellipse in the (R , N) plane 


+g = 1 
a* b 
dZ = - b* xX = - Tan L 
aX ae Z 
and 
Tan § = Z 
xX 


a: Qh0 = 
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thus a 
Tan§ = 2, tank 


But from the polar form of the equation of an ellipse 


= ab 
(a2 sin? S + b2 cos2$je 
or 
r cosé = a 
Se a rt 
(a2 Tan?§+1) 2 
be 
r sind = b 
(b2 cot “S+1) 2 
a2 
so that upon substitution 
*n 7 hy 
rp = a 
ee a ee eee 
(b2 sin? L + cos* L) 2 
ae 
+ b ith 
a ee 
(b2 sin? L + a“ cos~ L)2 


where: a = equatorial radius (Ro) 


o 
HF] 


= polar radius (Rp) 


Se 
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Computational logic 


=> ™> > 















Unit Vectors 


= Cosh eos RA. X + cosh sin R.A. g + Sink i 
al 


: - snRA RR + cas R.A. P ¢ Oo WN 
e -sink cosR.AR - gink Sin RAP ? cos LN 





Tracking Station Position | | 
a R a a 
8 (Sth) (uk + u,?) + (Fe +h) (uF) 
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2.3.1.4.2 Subroutine FQINGX 


Purpose: 


Deck Name: 


Calling Sequence: 


Input/Output 











ROTI 


Subroutines Required: 


Functions Requireds 


Approximate Deck 


Leneth: 


the true equator of d7*~ 
mean equator of 1950.0 (J.D. 74 33782.423) frame. 


CALL EOLNOX (TINE) 











FORTRAN 
Name 


TIME 


ROTATE 








computes the transformation matrices relating 
frame of reference to the 















Definition 


mean solar days since 
1950.0 (J.D. 2433282.42 





rotation matrix trans-~ 
forming vectors in the 
frame of 1950.0 to the 
frame of date 





rotational transforma- 
tion to convert vector 
in frame of date to 

frame of 1950.0 














nutation matrix 
relating vectors in 
the mean equator of 
date frame to those 
in the frame of date 
(used in GHA) 


MATHPY (matrix multiplication) 
ine) 


3 
cos 


= 
® 
—_ 


~~ rN 


1540 (Octal) 
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Description and Formulation: 


Over 2000 years ago it was discovered that the Verna’ Eauinox would 
move from east to west by 50", 2453 every year. This rn. .on is called pre- 
cession and is caused by the gravitational attraction of other celestial 
bodies acting on the equatorial bulge of the earth, If the earth were 
perfectly spherical and radially homogeneous, it would not experience any 
deviation from its mean equatorial pole. However, since the earth has an 
equatorial bulge, it experiences torques from the gravitational attraction 
of the sun and the moon, Due to the fact that the lunar orbital plane is 
approximately 5° oblique to the mean ecliptic, both the lunar and solar 
torques tend to align the equator with the ecliptic. The earth responds to 
this torque much like a spinning top responds to a torque. It precesses 
about the mean ecliptic pole. This precession is called luni-solar precession, 
Since the moon is so much closer to the earth than the sun, its contribution 
to luni-solar precession is approximately twice as much as that from the sun. 
The equatorial pole has an obliquity of about 23,5° so at the rate of pre- 
cession mentioned earlier, the equatorial pole would very nearly trace a 
right circular cone every 25,800 years. 


Just as the sun and moon cause the equatorial pole to precess, so do 
the planets of our solar system cause the ecliptic pole to precess; however, 
the magnitude of this planetary precession is very small and will be considered 
negligible in this discussion, 


"Total general precession" is the sum of planetary and luni-solar pre- 
cession and gives the changes in the mean vernal equinox of date from some 
epoch, Total general precession amounts to 50%, 2453/year and can be con- 
sidered uniform for practical use, This is the rate of westward rotation of 
the mean vernal equinox of date. 


As the equatorial pole precesses about the ecliptic pole, it also 
experiences further disturbances known as nutations, Free Eulerian Nutations 
are those which would occur if the earth were simply set in rotation and left 
to itself without any disturbing forces, Forced nutations are those which 
are caused by the changing positions in space of the sun, earth, and moon, 
which in turn cause variations in their respective gravitational attractions 
to the earth, 


The most significant nutation is the 19 Year Lunar Nutation. This 
nutation is causec by the precession of the moon's orbit, which is inclined 
about 5° oblique to the mean ecliptic, The line of nodes associated with 
these planes precesses with a period of about 18,6 years, The resuit is to 
change the direction or the small fluctuations in potential experienced 
by the earth-moon systems 


Other forced nutations include the Semi-annuai Solar Nutation and the 
Fortnightly Lunar Nutation., These phenomena are the result of the decreasing 
torque that the sun and moon apply to the earth as they approach the passing 
of the equatorial plane. Due to symmetry, the net torque, as one of these 
bodies passes through the equatorial plane, is zero, 
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The earth model that is generally used for the analvses of this motion 
is a rigid ellipsoid which is later simolified to an oblate spheroid. This 
model does not account for elasticity, fleidity and other physical properties 
of the earth; but it is sufficient to use for a fairly complete derivation 
of precession and nutation, It must thus be noted that the results of an 
analysis using such a simplified model are not exact, 


Improvement of the theoretical analysis resulting from the incorporation 
of observed data in the evaluation of the constants of integration and from 
the incorporation of more complete models of the earth in the analysis 
have, however, been effected. The results of these efforts based on formu- 
lations presented in the American Ephemeris and Nautical Almanac will be 
presented in the following paragraphs, 


Precession: 


Uniform precession is a rotation of the coordinate system defining 

the mean equator of date and may thus be represented by the matrix equation 
ci P P59 

where 7 0 denotes the position vector in the standard inertial reference 
frame (2ondamental plane ard principal direction are the mean equator of 
zero hours, January 0, 1950 and the corresponding vernal equinox) and 
where T_ is the position vector in the mean equator of date frame, The 
problem thus, becomes one of determining the elements of P, This step in 
turn is accomplished by selectiny 3 small parameters which relate the two 
frames, One such set is shown in the following sketch 


,mean Gales te of 


mean ecliptic of 
a 1950.0 





é ae ee ee 

va Ee - _ mé@an equator of 

Soo oy ho an 190.0 
Tn 


ie) 


ae one 
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The earth model that is generally used for the analvses of this motion 
is a rigid ellipsoid which is later simplified to an oblate spheroid. This 
model does not account for elasticity, fleidity and other physical properties 
of the earth; but it is sufficient to use for a fairly complete derivation 
of precession and nutation, It must thus be noted that the results of an 
analysis using such a simplified model are not exact, 


Improvement of the theoretical analysis resulting from the incorporation 
of observed data in the evaluation of the constants of integration and from 
the incorporation of more complete models of the earth in the analysis 
have, however, been effected. The results of these efforts based on formu- 
lations presented in the American Ephemeris and Nautical Almanac will be 
presented in the following paragraphs. 


Precession: 


Uniform precession is a rotation of the coordinate system defining 

the mean equator of date and may thus be represented by the matrix equation 
Fy = PBso 

where Yr 0 denotes the position vector in the standard inertial reference 
frame (Pendamental plane ard principal direction are the mean equator of 
zero hours, January 0, 1950 and the corresponding vernal equinox) and 
where fr is the position vector in the mean equator of date frame, The 
problem thus, becomes one of determining the elements of P, This step in 
turn is accomplished by selectiny 3 small parameters which relate the two 
frames. One such set is shown in the following sketch 
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Nutation: 


The relationship between the mean and true equator of date may be 
represented in terms of two small parameters as shown in the following 
sketch 


N 








P— Nean equator 


Yy €-€-S¢ 


t— True Equat or 


and the small parameters ( 6 wand S€ ) can be divided into long 
(AWand AE ) and short (dw and dé ) period contributions which 
can be computed as a function of a set of quantities defined in the 
Nautical Almanac. These quantities are given in time series as: 


OQ, = 129112790 -09052953922 p +020020795T +09002081T? — +0.000002T° 

@ = 649375452 4139176397 D —«- -9002131575T -0S00113015T* -020000019T% 
T’ = 208984399 +0°11190408 p -0°010334T = -0010343T* = -0°000012T° 

TP = 282908053 +02000047068 D +0°00G45S3T +0,0004575T2 +09900003T° 
L = 280908121 +0998564734 D +02000303 (T+T?) 


where D = Days since reference epoch (1950,0) (J.D. 2433282.423) 
T = Julian centuries past reference epoch 


Corresponding to these time series the Small parameters are: 


= 2eic 
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Foe 


ay 10° 


I Wx 104 


ae x 10+ 


dex 104 


Now since the 


ol 
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and since 


+ 


-(4798927 + ,0482T) sin. 
25800 sin 2M - 3.5361 sin 2L - .1378 sin (3L- T’) 
1 
60594 sin (L + J") + 0344 sin (ZL - (2) + .0125 sin (2T ~- 


23500 sin (L = 7) + ,.0125 sin (2L - 27") 


~95658 sin 2 €@ - ,0950 sin (2@ -f.) 

00725 sin (3@ - T') + ,0317 sin (@ +T*t) - 
0161 sin (@ -T! +) + .0158 sin (@ -~T' - 1) 

20144 sin (3@ +17! ~ 2L) - .0122 sin (3¢ -T' -1) : 
.1875 sin (@ -T') + .0078 sin (2¢@ - 2T') 

60414 sin (@ + T! - 2L) + ,0167 sin (2@ - 2L) 


00089 sin (4@ - 2L), 


2595844 cos ML - .2511 cos 222 
1.5336 cos 2L + ,0666 cos (3L - 7 ) 
20258 cos (L t+ T™) + ,0183 cos (2L - 1L) Ps 


.0067 cos (2T' - £1) 


02456 cos 2@ + .0508 cos (2¢ -dfh) 
.0369 cos (3€ -T') - .0139 cos (¢ +T') : 
0086 cos (@ - T'' + (2) + .0083 cos (g -T! - 21.) 


0061 cos (2g +T! - 2L) + ,0064 cos (3@ -T!' -N.) 


mean obliquity of date is given by 


2364457587 - .01309404T -,00000088T? + ,00000050T3 


E+ &€ 
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the rotational transformation relating the frames (X,, Yn? 24) and 


Br = Ty (-€) Tz (-Sy) Ty (E ) 7, 


eg a 
= N lm 
where T (  ) denotes a rotational transformation, the subscript denotes 


the axis of rotation and the term in the parenthesis denotes the angle of 
rotation, Expansion of this transformation yields N as: 


Ni) = cos éy 

Nyy = sin 5d Woos € 
N13 = - sin éy sin € 
No, =~ sin d Yeos € 


M, 


n22 = cor &Weosé cos€é + sin€ sin 


no3 = cosdWceoseé sin€ - sin€ cosé 


N31 = sin SYsin € 


N39 = cos d sin € cos€é ~ cosé sing 


N43 = cosdWsin € sin€ + cos€ cos é 


which upon substitution of E+ SE for € may be approximated as 
1 - Sycosé - Sysin e 
N =/5Weos € 1 -SE 
bWsin E d€ 1 


Combined Nutation and Precession: 


The true equator of data frame and that corresponding to the mean equator 
of 1950.0 can now be related by direct substitution of the results of the 
previous paragraphs. 


Yn = P Pso 
PT =N Tn 
bo a ~ > 
rp = NP Yeo = [ROTATE] rsq 
end 2 iis 6 
Poq = (NE)? , % [ROTINV] Pp 
- 3.53 as 
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Computational Logic: 
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Precession Matrix 


P= f (t) 


nutation parameter evaluation 
Oe Gers br =f Ct) 
Ayr dP, SE, dé = (0, ¢ vs Ly P) 


—om, 





Nutation matrix 


N=f (AW dp, A€, dé ) 











Total Nutation Precession Matrix 


ROTATE = NP 


RETURN 
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2.3.1.4.3 Subroutine GHA 


Purpose: GHA Computes the hour angle of the Greenwich meridian 
relative to the true vernal equinox of date in degrees. 


Deck Name: GHAN 
Calling Sequence: Call GHA (T, D, GH, DA, OMEGA) 


Input/Output: 


FORTRAN Math Common/ 
Name Name | Dimension | Argument Definition 


Universal time for the 
selected date at which 
GH is to be computed 





(sec.) 
ng D Deg 1 Arg Mean solar days elapsed 
sinceO Jan. 1, 1950 
(integral number 
0 GH Tr(t) if Arg Greenwich hour angle in 
degrees 
I DA da zB Arg Nutation correction to 
reference GH to true 
equinox of date 
6) OMEGA w ls Arg Rotational rate of the 
earth at the selected 
epoch 
Subroutines Required: None 
Functions Required: Sign 


Approximate Deck Length: 160 (octal) 
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Formulation: 


The hour angle of the Greenwich meridian relative to the mean vernal 
equinox of epoch T is given in the Nautical Almanac as 


P(t) = 1002 07554260 + OF 985647346d 
+ (2% 9015) x 10723a% + wt (mod 360) 


where d= the integral number of days past oh 1 January 1950 


t = time in seconds past oh of the epoch data 


sek CQLL7SO717 
1 + (5,21)10-34 


Thus, in order to be consistent with the measure of time in the remaining 
portion of the program (time reckoned from 1950.0)it was simply required that 
the whole and fractional number of days (Ds) reckoned from 1950.0 be changed 
by the difference 


d= Dso ~ -OTT 


and that a correction for nutation (da) be applied to compute the hour angie 
with respect to the true vernal equinox 


da = Speos é 
W(t) = Plt) + da 
~ 261 - 
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Computational Logic: 








w = ,OOKL78O7h2 
1+ 5.21 x 1074p 


GH = 100.0755h + .98564735DD 
+2.9015 x 10723DD* + wt 


GH » GH ~ X (36u)(Sign (1.,GH)) | 
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2.4 The Data Filter Group 


The reduction of the data provided by the preprocessor (PROCES) will 
be accomplished utilizing a minimum variance recursive filter first developed 
by Dr. R, E. Kalman (the reference and a description of the filter will be 
presented in KALMAN). However, because there are several distinct operations 
invoived in the process, it has been deemed desirable (from the standpoint 
of program modifications for other types of data, ease of development and 
checkout, etc.) to construct the filter in the form of several subroutines. 
This group of routines is the subject of the discussions which follow 
commencing with the driver (FILTER) and with the formation (KALMAN). These 
routines are discussed prior to all of the remaining routines, since they 
establish the mathematical framework in which the others can best be under- 
stood, and because they establish the notation to be employed. 


The interface of the differential corrections program and its pre- 
processon is also embedded in this group. This interface exists in the form 
of a routine (DATAPE) designed to read the specially prepared data tape 
and load the smoothed observation vector and identifying information into 
memory. The operation of this routine has been checked (as has the operation 
of all other routines in the program) to assure compatibility between the 
two programs. 
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2.4.1 Subroutine FILTER 


Purpose: FILTER is designed as the driver for all routines in 
the Data Filter Group. It also serves the function of 
computing the observation vector (observed minus computed 
residuals). 


Deck Name: FILT 
Calling Sequence: CALL FILTER (NUMB, R, RD, A, E) 
Input/Output: 


FORTRAN Math Comnon/ 
7/0 | Name Name Dimension | Argument | Description 
I NUMB - 1 ARG 


Number of the station being 
checked at this time in the 
routine TRAK 


I R - 1 ARG The computed values of 
RD Va Hi ARG range, range rate, azimuth 
A A A. ARG and elevation relative to 
E E 1 ARG the tracking station (based 
on the optimum estimate of 
the trajectory at this time 
(Km, Km/sec, rad, rad) 
IT ODATA = 3 WRK(65) Observed values of range, 
range rate azimuth and/or 
elevation (no more than 3 
pieces of this information 
at a time) (Km, Km/sec, rad, 
rad) 
I ITYPé ~ ak WRK (64) The fixed point index which | 
identifies which of the six 
possible types of data is 
being processed 
; , ? 
Subroutines Required: MEASURE (computes 25 ) 
ERROR (computes @ ) 
KALMAN {computes ¥ and Pp) 
UPSTAT (updates problem) 
DATAPE (provides next data point) 
MAIN (deck name of main program) 
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Functions Required: None 


Approximate Deck 
Length: 205 (octal) 


Description: 


Subroutine FILTER serves as the driver routine for the complete filter 
package. It is designed to funcvion in such a manner that all of the informa- 
tion required by KALMAN is available when the call is made. To be specific; 


1) The firs: step made is the identification of the type of data being 
processed and the construction of the ordered set of observed minus 
computed residuals. 

“= T 

2) The next step is the computation of the matrix uM? CopstT = =f ) 
in subroutine MEASUR and the definition of the weighting matrix in 
subroutine ERROR 


3) At this point KALMAN is called and the state vector (X) and the 
covariance matrix for the error in X are estimated. 


4) UPSTAT is then entered to determine if the trajectory is known 
to a degree sufficient to allow-up-dating. If so, new conic 
elements are computed and written. 


Upon return from UPSTAT the problem has been reduced by one of the smoothed 

data points (7 words — TW, TF, TYPE, STATION, OBSERVATION VECTOR (3) ) and 
DATAPE is entered to determine if additional data are available. If so, the 
the next data point is read into memory and control is returned to the program. 
If there is no additional data, exit from the program is accomplished by calling 
the $IBFTC name of the main program (MAIN). This sequence is utilized at this 
time since it allows other cases to be run in a sequential manner, however, it 
is noted that a modification of this procedure may be required for other 
systems, 
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2.4.1.1 Subroutine KALMAN 


Purpose: 


Deck Name: 


Calling Sequence: 


Input/Output : 


| FORTRAN 
I/O | Name 
I 


OBVEC 


I M 

I PHI 

I OBST 
I Q 

1/0 | STATE 
1/o | P 


To obtain the minimum variance estimate of the position 


and velocity relative to 


an estimated trajectory and to 


produce the covariance matrix for errors in the estimate. 


KALM 


Call KALMAN (OBVEC, M) 


Math Conmon/ 
Name | Dimension Argument 
Y M ARG 


2 Ll ARG 
Pltt)| 6x6 STT(1) 
mit) | 6X3 STT(37) 
Q(t) | 3X3 STT(55) 
X (#) 6 STT(64) 
| 
P(t) | 6X6 STT(70) 
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Description 


The vector composed of 
observed minus computed 
residuals (from FILTER) 


The dimension of the observa 
tion vector 


Transition matrix relating 
errors at two successive 
data points (from STMAT) 


Matrix of partials of obser- 
vations with respect to the 
state (from MEASUR) 


Covariance matrix for contri- 
bution of errors in observa~ 
tions and of errors in 

station location (from ERROR) 


ar 
ee vector § ay f where 
r = 


r- To 
— =” ole 
sub o = reference 
Covariance matrix for errors 
in the estimated state 
vector 





SID 65 1203-1 


Subroutines required: TRANSP (matrix transpose) 


MATMPY (matrix multiplication) 
MTINY (matrix inverse) 
SUBMAT (matrix subtraction) 
ADDHAT (matrix addition) 
Functions required: Wone 
Approximate Deck 
Length: 1056 (octal) 


Formulation: 


The theory of Kalman's recursive minimum variance data filter was first 
oresented in a series of papers (e.g. Ref. 1) in which the author jeveloped 
rigorousiy the form of the resultant estimate by employing an orthoconal 
projection lemma derived as a portion of the naner. Since this computational 
algorithn is utilized in the digital progran being discussed, its form must 
be discussed. However, since complete descriptions of the development are 
recorded (Ref. 1), only a summary of the steps required will be presented. 


The basic assumption of this nrocedure is that the optimal estimate of the 
state vector (X) for the system (in this case the vector oe ) is of the 
forn. q* = M8 Qi Yi where Yi. for this discussion are the com- 
penents of thé” ‘observed r:inus computed residuals (Y = MK), where the systems 
equations are 


x(t)= (6,4) x(4) + u(£), 


where y (€,4& ) is an mbyx matrix of time-dependent coefficients and where 
u (t) is an independent Gaussian process (u will not be utilized 

in the estimation procedures). Lest this assumption be questioned on the 
grouids that it is excessively restrictive, note is made of the fact that Kal- 
man proved that tne results obtained with this model can, in general, be 
iiaproved a with non-linear estimation only if the errors in the data and/or 

in the state vector for the system are non-Caussian. Further, to achieve the 
improvezent, at least third-order distribution functions for the errors are 
required. 


Now, attention is turned to the problem of obtaining the optimum estimate 


of the state (x* (e) ) at some tine t by utilizing only the last estimate of 


‘Une state (¥* (4-1) ) and the observation at time t. This statement of the 


roblen leads es ne Naeue for. of the computation algorithi 
Pp ‘ 8 


RE)? le es) R*(E-de K(E)[ YE) 
-m(E) gy (6&1) ¥* (4-y)] 
=p (6, él) X*(tl) + K(t) ag (f) 
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" where the first term is the estimate of X(t) obtained utilizing data a. quired 
prior to the last data point but propagated to time t. The second term is 
a weighted correction which is determined by the difference in the observed 
and predicted values of the observed minus computed residuals. 

The task is now to define the optimum linear gain K(t) (optimum here will 


be taken to mean in the sense of minimum variance) to be utilized in the 
algorithm. This task iill be accomplished by adopting the notation 


y(t) = M(é)x(é) +€ (6) 


x(é) = ¥* (4) + WCE) 


where € (t) and 7 (t) are Gaussian errors and where the notation  (t) means 
the error at t based on all data processed prior to t. Expanding these 
identities 


y(t) = M(t) P(E, é-1) R*(E-/) 
+M (é) P (4 4-1) He) + E(E) 
014) = [p (6, 4-1) ~K(E)M( ¥(6,E-D) (E) 


- K(é)e(é) 
and noting that the optimum estimates for X satisfies the conditions that the 
covariance matrix forg@ (t) 4 7(t) is minimized, or from the referenced lemma, 
that 
E(4¥ (4) ¥"Ce)) =o = Expected value of ( ) 


ora upon substitution and expansion (after the independence of ¥ and 
and of E and X are assumed). 


[y(e en) - K(e)M (8) o(&E-“)) Ely (é-)7 (eee e-1)m 7(é) 
- K(k) & [e (6 €"(e)] 0 


At this point in the development, a slight change 1 the notation is 
adopted in that 


P(é-l) = Eg (e-) ”(é-/)] 
Q(é) = Efee' CA 
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and it is noted that if P (t-1) is positive definite then the product 


M(t)¢ (4,41) PU-l) G7 (é, EME) 


Will also be positive definite provided that the observables (components of 
y) are linearly independent. Thus, the product is invertable; and the 
optimum weighting matrix is 


K(t)= Blt, t-1) P(t-!) p(t, t-1)M"(t)- 
[ce]? (88-1) PCE p"(E, ey) MME) #Q (EN) 
= P(t) M™(E)[ M(t) PLE) MU) + Q(t) 


Now since the distributions of the errors in the estimates can be prescribed 
at time zero and since the covariance matrix for the contribution of the 
stations to the uncertainty can be defined at any given time, the first 
weighting matrix can be computed. However, to reduce additional data points, 
the relationship defining the covariance matrix for the errors in the estimate 
after the reduction of the data point (P (t) ) must be defined. This 
matrix may in turn be evaluated from the definition of P (t) 


pP*(t) =E/* (t) r*"(t)) 
= p*(t,t-1) O(t-!) p*7(t, £-1) + KL) QU) KL) 


where P (t,t) = Plt; t-1)-K(t) M(t) p(t, é -/) 
Filter Description (operational characteristics) 


The minimum variance formulation presented here is a simplification of 
the maximum likelihcod estimator in that it assumes the errors are normally 
distributed. This, however, is the only valid objection voiced to the use of 
this procedure since in its general form the formulation includes most other 
data filters. (For example; weighted least squares - the simple case where 
the components of y are uncorrelated). The recursive nature of the filter also 
provides a distinct advantage since non-linear effects resulting from errors 
in the equations of motion and from non-precise relation of the errors at 
various points along the trajectory are minimized by the fact that time 
required for data accumulation is itself minimized. Further, if the trajectory 
is updated at each of the data points (assuming that the elements of P (t) are 
sufficiently small to make this practice feasible) the linear system can pre- 
dict the behavior of the non-linear system to a very good degree. Finally, 
due to the recursive nature of the filter, the order of complexity in reducing 
any number of data points is constant; and problems of loss of numerical 
significance arising from maxipulating large arrays of numbers are drastically 
reduced, 


Ref. i Kalman, R. E., "A New Approach to Linear Filtering and Prediction 
Problems" Journal of Basic Engineering (March 1960.) pages 35-45 
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Computational Logic: 


’ KALMAN 










update covariance matrix of Ar,Av 


Py py! 


construct weighting matrix 


w= Put [opt + Q] . 





weighted error vector 





OB = OB - MY fstate? 



















update state vector 


state = W.0B + 9 {state} 







atone 


update covariance matrix of Ar , &V 


P= Pp —- WMP 


RETURN 


(ao 
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2.4.1.2 Subroutine STMAT (State Transition Matrix) 


i med Purpose: 


Deck Name: 
Calling Sequence: 


Input/Output: 


FORTRAN Math 
Name Name 
F 





to provide the 6 X 6 matrix of partial derivatives o 
the state vector at some arbitrary time with respect 
to the state vector at an earlier epoch. 

STM 


Call STMAT (TIME) 


radius vector in frame 
of date at time t7 


] 


velocity vector in frame 
of date at time tj} 


radius at to 


velocity at to 


Subroutines required: TRANS (conic Transition matrix) 


INVAO (analytic transition inverse) 


Functions required: None 


Approximate Deck 
Length : 


310 (octal) 
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Description: 


The trajectories for the vehicles of intcrest to this stuiv are nearly 
conic. Thus, partial ierivatives evaluatei for the nowinal conic trajectory 
will agree well with those obtained frou. the truc trajectory by nore 
elaborate means such as the integration of the adjoint eauations. Yor this 
reason, the conic representation will be utilized to construct the c X 6 
matrix of partial derivatives though two slight uodifications will be 
employed to assure that the conic and truce matrices agree as well as pussible. 
The first modification is that the position and velocity vectors usel tu Je- 
fine the partials will correspond to those of the true trajectory at the tine 
from which the errors are propagating in this case, (the last data point) 
rather than the corresponding point on the conic references trajectory. The 
second modification is that the true position and velocit; vectors at the 
tine of the present data point will be utilized to obtain the inverse of a 
second estimate of this matrix by solving the conic problem backward in tine, 
This second matrix will then be analytically inverted, utilizing a special 
propesty of this matrix (developei in INVAO) to obtain the desirei partials. 
These two matrices of partials will not be identical because the conic 
trajectories utilized to define them differed in regard to each of the six 
constants of integration (the result of oblateness, drag --- forces). 
Further, they will differ from the true matrix. However, a first-order esti- 
mate of the effects of these forces can be included and the agreement with 
the true solution improved by averaring the two separate solutions for the 
partials, 


Computational Logic: 






Y= Y(rtran, vtran, t) 






~t 
O = (rvec, vvec, tt) 


3 
i 


US teen 
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2.4.1.2.1 Subroutine TRANS 


Purpose} TRANS computes the matrix of partial derivatives 
of the state vector at an arbitrary point on a 
conic trajectory with respect to the state 
vector at another epoch. The variables utilized 
in thig analysis are well defined for all conic 


motion, 
Deck Name: TRAN 
Calling Sequence: Call TRANS (R, V, T, PHI) 


Input /Output: 


Definition 


radius vector in cartesian 
coordinates at T = o (Km) 


velocity vector in cartesian 
coordinates at T = 0 (Km/sec) 


time at which partials are 
desired (relative to point 1) 


aX2 


matrix ae 


gravitational constant 


Subroutines required: SEARCH (solve analog of Kepler's equation for 
position as a function of time) 
Functions required: AMAG (vector magnitude) 
Dor (dot product) 
SQRT (square root) 
COs 


SIN 

COSH (hyperbolic sine) 
SINH (hyperbolic cosine) 
DERAQ 6 = 0 6% 

= 1 i =Y ) 
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Approximate Deck Length: 820 (decimal) 
Formulation: 

The equations of conic motion in terms of the -o-cclled universal 
variables of Dr. S. Herrick will be utilized to develop the partial deriva- 
tives of the components of position and velocity at any given time (on the 
conic section) relative to the components of position and velocity at some 
other arbitrary epoch. This task will be performed utilizing a formulation 
valid for a non-rotating coordinzte system and will be based on the develop- 


ment presented in the discussion of the reference trajectory. The required 
expressions for this analysis are: 


PPE +98 

3-fF + QS, 
where f= / - Sp, 

g*T-U 

fe -S/re 

g=/-C/r 


= inertial position vector =X z* + yy +zz 


4h 


_ ~ AX . “oA 
= normalized velocity vector = V/VZ = XX" + yVrzz 


uy 


A 
xX - is the X unit vector. This notation is adop*nd to avoid confusion 
with a variable to be defined subsequently. 


€=a(/-coox) 


4 . 
G=a*(x-am xz) 
3 =a* (mx) 
=- 5K 


X= £-E, ¢ elliptic motion ) 


~22Q.. 


STM GE 12021 


ae panes a4 a : = ee OT ae, ET iy OER EYRE 8 SERRA ree | 
a llreaeacecea ana ae 


a na rata ' mS ee = . Z 


pea ae (hyperbolic motion) 


The first step in obtaining the desired partial derivatives invcelves 
differentiation of the equations for ¥ and 6 with respect to the components 
of To and Bo: This task will be drastically simplified if full advantage 
is taken of the similar form of these derivetives at the outset. Thus, a 
shorthand notation will be adopted in that u and v (u and v) can assume 
the values of s, y, and z (x, y, and 2) 





2, =f5,, +U, Fe 
“eg nee 
a FO uy UF + Use 
ee ar a U 


where 6y,°O &#v 


=/ uszy 


Thus, the problem has reduced | itself to one of obtaining the derivatives 
of f, &) f and & with respect to r and 8. This task will in turn be 
simplified if a set of intermediate parameters is selected, since the x...8 
do not appear explicitly in the equations for f...g. The set to be utilized 
is suggested by the equation for the magnitude of r in this set of variables. 


rs +D, S+(1ty a)Ee 
=F (5,D, >) 


Having selected the intermediate variables the next task is the 
differentiation of f, g, f and &e 





For f 2 re 
Oi, Cas = 
an, nt 
Ft we 
ID, r, 9D, 
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ON sc: doy aC 
Ga i. ae 
For ¢ x 
ag . ar_ 2d 
04, 9% ah 
dg . 97 _2U 
3 3D, 20 
ag . ar_ al 
as a 
For f a 
> os ) 
= —_—_— tr 
oF (rre) oe +S(r °F, 
dr, (Fi) 
oF _ -(rh)znt Se Sp. 
a, Cr te)* 
. és a or 
ao (rrm)* 
and for g a 
. ec A or 
4 = "oR tS Op 
aN re 
; gE a 2 
% : 50. as 
2D, r& 
20 ac a or 
@  ~tan to Be 
ax r2 


ea AA 
Now attention turns to the derivatives of C, S, U, T etc. with 
respect to ry, Do, ®. 


for G 


for 0 





Q 
OC» 


> 'O> 
ce) L ot 


% +12 Q|o Y 


ey 8, 
tt 


NY 





. aX a 9% 
adn X On Ss an 
a Or 
g =e. 
aD, ' 
3 ga HOM 5G = S gor to 
az 


i) 
ot 


ax 
a® ( 1-co0 x) rig 


Oa 
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90 2a ax 
90, © a 0, 
aU a ag ty OQ 
ga So Get BlX-ainxja* Fe 
A aX 3a 
=C aqtgZue 
for § 
aS = may 2 
; ol, 
af. ax 
30 Coo XL 3D, 
as az 1 5 ¥ 
—— — a 
a COO X Bot +B One 
for v 
I? _o 
a4 
oT 
2" :0 
Do 
¥ 
a 2B 
an? for r 
or os C gc A 
co Hs pulled —- +X 
ar, / D ar * 3 oF +aXC 
er A 23 r ac 
oo” $+#D, 2D, tGy aD, 
ar as ac Bf. 


The final set of partials required is now recognized to be that of 
with respect to r,, and Do, and o& . This set requires the eguation for 


Vyas, 


time (analogous to Kepler's equation) be differentiated as follows: 


Ja CET eee Orc $e) 0 





ax 
for = a a “A 
on A ox OC “~ ou 
= —— D -OU rll t+ Xr, ) — 
OrX tle FP 2 on? (/ 0) ar 
“Rea Ser S 
& ae 6. 4 on 
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me = + - BRT TT peer errant tl Re . aaa SS 


9% __§ 
an r 
ax on ac goa 
On R>z- +Cr 3B 
for 0 aD, +D; ye Co FG, 
70, or 
ax Z ns 
for gy J si Oe ac a 4 
Get a oi haat » 5a t So 50, + U% 
= Ol+ %0.a)+2,(Ca)+ 2% zr) 
OX ok ae 
oe 2-11 (ca +6) +F(OE-1 x) 2 - 


Now, substituting back into the previous exoressions and collecting 
terms, the derivatives required to compute the Par beee. Of 185 Fs. ey 
with respect to the intermediate parameters are: 
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ay 71 t Oo(-Ceox )tG(-p)raé =2F 
ss Se). 4 
a7 8 Al cou x) #Cy (- : Jezt¢ 
er a 
BO * Di. Set C, Cut % = Ve 


Finally, the required partials of f, g, f, and &, with respect to 
Yo, Do, and ox are: 


Fer F Lee Mele 
% ne rr 
af =e 
9Q rr 
of ! 
aa *r, Sm 


For q 


See ¢ - a 


ba FR Oe rer re 
Garg aie dg 
3%" irs +Crnf] 
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The only remaining steps at this point are thus to provide the deriva-~ 
tives of the intermediate set of varameters with respect to the components 
of Lo» and Bos to construct the derivatives of f, g, f and & with respect 
to To and &, o» and to relate the results to the parameters Y and Vv. The 
first step ¢s accomplished by referring to the definitions of ro? Do, and 











3 3 
R 2 a0 ~ R 
‘0 -2 a3 70 tes % 
an X& 2h 26 
ay, % ? 9S; 
D,=& «8, 
aD, ID _ 
ax, 3é 3 aS; ei 
é 
lee —_. z 
a= S$ “Sop, 
Oe 5A OE Gs 
9X; “4 7 aS; ¢ 


The second step is accomplished through the medium of the chain rule, i.e., 
of oF a af ae af ga 


QZ, Or a4; * aD, ax; ” 2a AX; 
etc. 


and the third and final step is employed to remove the normalization factor 
applied to the components of the velocity. Since 


ate 











S ie 
ds = Ba ; 
the desired matrix is 
ar | Jt OF | 
dr LN ey ae, tem dr, 
dv oss as dv, 
ol a%~C«Sd as, 
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Computational Logic: 


Subroutine 
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Call Search 
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Fllintic 
Mation 


G=/(-/m)* (x-sinx) 


Mation 
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f= /-/r, 
ae aa 
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U* = - 4a (2/6 -0) 
ot =~ Ya (2% -€) 
Uan= 2 (ic*- 3u*) 
Ca = 2(20-2c*) 
Xe = fo Ut Do Ca + Co Va! 


Kz ]eac 
SS Ure Ux 
Fa ® 10 tO. Sa tC, Cx 
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2.4.1.2.2 Subroutine INVAO 


Purpose: to produce the analytic inverse of the 6 X 6 matrix of 
partial derivatives referred to as the state transition 
matrix (i.e., the matrix zh where ail vectors are 
expressed in cartesian coordinates) 


Deck Name; INVA 
Calling Sequence: Call INVAO (AO, AOTZ) 
Input/Output : 


FORTRAN 
1/0 Name 
I AO 
Oo | AOI 


Subroutine required: None 


Definition 


the matrix of partial 
derivatives of the 
state vector (x) at 
"t" with respect to 
the state vector at 
"to" 






Math Common/ 
Name Dimension | Argument 
6 ARG 




















the analytically in- 
verted matrix 









(en 


Functions reouired: None 
Approximate Deck 

Length: 160 (octal) 
Formulation: 


Consider a linear system (expressed in cartesian coordinates) described 
by the following equation 


k= Alt) X(t) +B F () (1) 
X(t) 3 even-ordered state vect d of a set 
Variables and thelr derivatives "'* Vector composed of a set of output 
~306- 
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we 

F(t) is a forcing function ( = 0 for the analysis to follow) 

A(t) is a coefficient array for the system composed of square, siTwtsic; 
even-ordered subarrays of the followins form 





° 1 A, 
A(é) = Sire ae (2) 
Aas fo) 
het 22 “|r - at" 
OF r are 
B(t) is an array relating the sensitivity of the system to the forern, 
function - 
and the solution for F =o 
— ara " Zaz —s 
X(Q)= Y(é4) X(4) = BD \M by) 
1 


If equations 2 and 3 are substituted back into equation 1, the followins 
identity results 


g Ge oO L E Z| 
Yar B2\" ‘ oy Y, Yoo (4) 


and equation 4 yields upon: expansion 


%* %, (5a) 
¥, > (5b) 
g = a ZY, (5c) 
he = 8 We Ca) 


Zovuations 5 may now he onerated on to produce an equivalent set of 
differential equations. This operation is performed as foilows: 


r 


ee 


27 Sgr ae | a got We: 
tar th th 32 Se * BSB 


E80 7= 
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cee ne 


or 


2 (4 Rh BH) 0 

CA ana i ak (6a) 
sinilarly 

Pao" 0 — bq" $y, 7 Ca (6b) 

Pa Ga “%, 2,76; (éc) 

G," $. ~ $y" $= (64) 


Finally, the results of the intesration can be restate! in matrix 
notation as 


= 
a ° ee % t Ca C, 

r . 7 
“Zao” $e, a, -C, Gy (7) 


and the constant arrays resulting from these intesrations rar now be 
evaluated by substituting the initial conlitions 


GD eg LO). ee 
%, (o) = a (0) +O 


This step produces 


C¢,+0 
Cx Sh, 
Coe 0 
Cyr t 


~208- 
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ani reduces equation 7 to the identity matrix. But, the only naetrix which 
ean be utilized to reduce an arbitrary scuare matrix to I is its inverse. Thus 


4 Gia 2G 
AEs ap : (8) 
) A) oF 


fquation & is inaportant for zeneral linear systens in that it provides 
an analytic means of constructing the inverse transition matrix directiy fro 
the elenents of the knowm transition matrix by rearrange.ent of terus ani the 
change of a few signs. It will be utilized in this progra.. in conjunction 
with a correction to the basic conic transition matrix Jescirbeld in STAT. 

In conclusion, it is noted that the true meaning of the tezms Ayo and 
A21 (equation 2) was never employed, and only sjmmetry is required. Thus, 
there 1s an imnediate seneralization providing that an arbitrary system. can be 
described by equation 2. It is also noted that another approach to the 
aerivation is suSgested in Ref. 1, along with a discussion of several relatel 
problems. 
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2.4.1.3 Subroutine MEASUR 


Purpose: To construct the matrix of partial derivatives of 
the observables with respect to the state vector 
(i.e., M= ot ) where the observables may be 


1) range 

2) range rate 

3) azimuth and elevation 

4) yrange and range rate 

5) range, aximuth and elevation 

6) range rate, aximuth and elevation 





Deck Name: OBSN 

Calling Sequence: CALL MEASUR (M) 

Input/Output : 

| FORTRAN] Math Common/ 

1/0 Name Name | Dimension j|Argument Description 

I ITYPE = a WRK( 64.) Fixed point number which 
identifies the type of data 
being processed and hence, 
the type of matrix desired. 

3 WRK(111) Radius vector in frame of 

I | RVEC F date 

I | VWEC ¥ 3 WRK(114)  |Velocity vector in frame of 
date 

IT RTRAK Ty 3 WRK(95) Radius vector of tracking 


station in frame of date 


I X (Y,Z) |U (E,N) | 3 (3,3) WRK(108) Up (east, north) unit vectors 
WRK (1.24 ) at the tracking station 
io 


AKG Fixed point integer which 

| identifies the number of 
pieces of data in the obser- 
vation vector (1, 2 or 3) 





wat on 


aT) 6 3202-1 


FORTRAN Corsaon/ : 
Name 3 Areruzent Oe ad ohh ees 


isfeae es rad i ° 7 
aVray 2? yaxtcal 
vatives cf the obser- 


vations wit, 
to the state 


OTGA 





Suoroutines required: CLRO33 (cross product) 

Tunctions required: AMAG (vector macnitude) 
er (Jot product) 
SCRI (square root) 


Approxirate Deck 
Length: 350 (decimal) 


Formulation: 


Partials for range, range rate, aginuth and elevation with respect to 
the state vector can be obtained from the followine equations: 
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where: 


b=) 
1 


= rango 


' 
f; 


= range rate 


€ 


Ag = azimuth 


DL = clevation 
nN A on 
U,E,! = up, east, north unit vectors 
~ = vehicles position in equatorial frame of date 
tn = nominal position on reference orbit 
Py; = tracking stations position vector 


when it is noted that for the purpose of differentiation, the nominal posi- 
tion vector and the tracking station position vector are constant, i.ec., 


This set of operations has been performed, and the results of the 
analysis are presented below: 


1) Partials of range 
2 
R°=XP + Ye + 2,7 


SE (GB + F)-(B Boe 


2) Partials of ne 


Zz ,0,0,0) 


3) Partials of Azimuth 
8*2 (7-2) + (F-R)* 


RP) NER) 


st $v 


aAg = (2Ag > tha) = [ &( 
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4) Partials of Rlevation 
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Computational Logic: 
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2.4.1.4 Subroutine ERROR 





e Purpose: To compute the weighting matrix (Q) for the data point 
p (one of 6 types) from uncertainty data provided per- 
taining to station locations errors and recording 
instrument errors 


Deck Name: EROR 
Calling Sequence: CALL ERROR (NUMB, M) 
Input/Output: 
FORTRAN | Math Conmon/ 
I/o | Name Name Dimension } Argument Description 





1 ARG Number of the station ~- this 


number corresponds to the 
sequence established in 
INPUT (comes from TRAK) 


1 ARG Number of pieces of informa- 
tion included in the obser- 
vation vector (from MEASUR) 


1 CON(1) Earth's equatorial radius 
(Km) 

1 CON(2) Earth's polar radius (Km) 

uf CON(7) Earth's spin rate (rad/sec) 

3 SDA(51) Variance and latitude, 


longitude and altitude for 
the station 


M SDA(141) Variances on the observed 
data for the equipment 
utilized 

6X3 STT(37) |Output of MEASUR 

3X3 STT(59) Combined weighting matrix 

a | WRK(64) Type of data being processed 
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Common/ 
Argunent Description 


WRK(105) | Geodetic latitude of the 
station (rad) 


WRK(106) | Iongitude of the station 


(rad) 


WRK(107) | Altitude relative to 
reference ellipsoid (Km) 


: ;) WRK(108) |p (East, North) unit vector 
at station 





. Subroutines required: MATMPY fates multiplication) 

TRANSP matrix transpose) 

Functions required: SIN 
Cos 
SQRT 

Approximate Deck 

Length: 627 (octal) 

Formulation: 


In order to construct the matrix of uncertainty required by the Kalman 
filter resulting from station position errors and errors in the observations, 
two assumptions are necessary. The first is that the groups of errors from 
the two sources under analysis are statistically independent, and the second 
is that the equations desc.sibing the processes are linear. Under these 
assumptions the resultant uncertainty matrix (Q) is 


Q:=@Q,+ Q, 
and attention can be directed to each source independently. 
In order to define Q], it is first necessary to consider the manner 


in which the specified set of errors affect the problem (consistent with the 
assumptions previously outlined). This relationship is 


4b 
AA 
Ad# 


Q(08s) = sees) 
oL,d,H 
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oe a ae - Ree re ee peep <I an AE oe - Dyn eerennemerneecmee enemies eM met See Se wide, 
ERC or ROE oo : + 


es 


_ 


aa 
cere TE 


= [COEFF] ax} 
AH 


noha oe cans 


which leads directly to 


a 


E (A (0as) A” (0es))= |courF] F (f an fab A 4H?) [coErF| 


Q,= [cosrr] g- [cozrr] 
where @ is specified. 


The problem to which this analysis must address itself is thus the 
construction of the matrix [COEFF] . This development is simplified 


considerably, since the chain rule can be employed to produce the desired 
results as follows: 


come = [ 2 (638) || 7 », | . M33 a 


=[0as}[ DERIV] 


wher [oBs] T is the output of subroutine MEASUR and where X is the state 
vector. For this reason, subsequent analyses will be directed to the deriva- 
tion of the matrix DERIV. 


The first steps in relating uncertainties in the position and velocity 
of the tracking equipment and the corresponding uncertainty in the satellite's 
position and velocity are the definition of the observed quantities 


pe FHF 
is we wk 
V,= V- Vr 3 


where: the subscripts r and t denote relative and tracking station, 
respectively, and the description of a coordinate system in which the vec- 
tors will be defined (in this case the true equator of date frame). 


7p * (Rey H) Scook coo 8 + 600 Lowa § 
CG 

* (fe oH) Jeouel zZf 
V7 W(2 xk) 
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ee wee Cee aie Oe aa Oe ila 7 ii ae ES tel 
Op Sp ror penser DEEL Bars call ee : = 








¥. =W(K,. | aceasta Kr cool cooA P| 
or 


=-w4 Rewx,YrOZ 


where Re 


3 


equatorial radius of the oblate spheroid 

polar radius of the oblate spheroid 

a function of latitude 

cos® lL +(K+/Re)* sin? L 

geodetic latitude of the station 

right ascension of the station 

longitude relative to Greenwich ( A ) nlus sidercal time 
of the Greenwich meridian 

spin rate of the earth 

unit vector toward the vernal eauinox 

unit vector along the spin vector of the earth 


re) 
toe 


owl 
f 


rape 
t i i 


Now if the vehicle's position and velocity are held constant for the 
purposes of differentiation, the sensitivity of the state vector to errors 
in the tracking stations position and velocity can be obtained by computing 
the partials of the relative Y and ¥ with respect to errors in L, A, and H 


ar = -dr, u=LA,H 
ou du 

Y ats Vy as LLA,H 
ou ou 


However, before these derivatives are evaluated, it is noted that since 
right ascension can be expressed as 


A> W(é-6)+#d 


and since sidereal time is uncertain to a degree much less than that associavel 
with A , dA Ss dad. 


Now 

OX, --(Re 4H) tien Cooh+ lool cooA(-Re\ dc 
aya (Be +H) (4) DL 
ok = “(Ke +H) CooL ain A 
oh C 
dX - Cook cooA 
dH 
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~—e 


where 


and 


ae. -(Re +H) aun ain + cool en A(R ) de 

dL ape e 
r =/kKe A 

at = (Be + H) cook coo 








C Ke Rec? L 
Zr. O 
dX 
0Zr. amb 
dH 
£86 -wxmk wok [/- (Re\* 
c* GL 2? [ (@) J 


FUNCT (in program language) 





aXr -W a, 
dk OL 
OX, = -WdYs 
gy doa 
aX, = - Wat, 
dH aH 
av, = ¢ Gt dX, 
dk OL 
Chee +@ aXe 
dA o> 
OP, =4W dK, 
oH, gH 
0270 “ek, 
ou 


With the evaluation of DERIV, the array COEFF is known and Q] can 
be computed once data for the errors in the variables L, A , and H are 


specified. 


Attention thus turns to the second source of errors to be con-. 
sidered; namiclzr, those resulting fron observational uncertainties. 


Assuming 


that the observation vector for this analysis is ordered in the same fashion 


TT 
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as for the previous analysis, Qo is by definition 
Q, = E(ee’) 


E = Yreue - Youserven 
and is provided as input data for each station. 


For the present, both @ and Qo have been assumed to be diagonal (i.e., 
the errors involved are uncorrelated); and provision has been made in sub- 

. routine INPUT to store only the variances. However, should station calibra- 
tion and survey indicate that a significant degree of correlation exists 
between the elements internal to these arrays, other terms must be added. 
This step requires only the modification of INPUT and ERROR. 


Computational Logic: 







Subroutine ERROR 


T= STATION ERRORS 





[Dea] 2 dF V 
dL dH 


| (Obs) -C -[M)} x[ DER!V] 


OL,d,H 
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2.4.1.5 Subroutine UPSTAT 


Purpose: 


Deck Name: 


Calling Sequence: 


Input/Output : 
FORTRAN 
I/O | Name 
1/o |STATE 
I P 
I ROTINV 
0 RCONIC 
0 TCONW 
TCONF 
0 RTRANS 
VTRANS 
o | PTRANW, 
TTRANF 


Math 
Name 


X (8) 


P(t) 


> a 


Tose 





To determine if the uncertainties in the estimate of ; 
the state of the system are sufficiently small to Basal 
justify updating the position (velocity) by adding 

the state vector to the position and velocity vectors 

on the nominal trajectory and to write the results of 

the data reduction problem. 


UPST 


CALL UPSTAT 


Common/ 
Dimension | Argument Description 


The state vector for the 
system (Km, Km/sec.) 
6X6 Covariance matrix for estima- 
tion errors in X (+) 

i ie Sa Matrix relating coordinate 
frame of date to frame of 
1950.0 


Position and velocity vectors 
for new conic reference 
trajectory (Km, Km/sec.) 


Time in days at which recti- 
fication of reference conic 
occurred (days relative to 
1950.0) 


Position and velocity vectors 
in frame of date on true tra- 
jectory at data point just 
reduced (Km, Km/sec.) 


Time from which errors will 
propagate for next pass 
through filter (days relative 
to 1950.0) 


1,1 





A227 
Ss™ 6& 1202-] 


ae! 


were 









Description 


Whole and fractional part of 
the present day relative to 
tk epoch of 1950.0 (J.D. 
21,33282.4,23) 


| Common/ 
Argument 
WRK(50) 
WRK(51) 


Dimension 





















| WRK(111) 


Position and velocity vectors 
WRK(11,) 


in the frame of date 
(Km, Km/sec.) 
















CON(1) 
CON(6) 


Fquatorial radius and gravi- 
tational constant for the 
| Farth (Km, Km3/sec*) 











CON( LL) Output tape nwnber 


Subroutines required: MATMPY (matrix multiplication) 
ELEMEN (computes conic elements) 

Functions required: DOT {dot product) 

Approximate Deck 70h (octal) 

Length: 

Discussion: 


UPSTAT (Update State Vector) is designed to determine whether the un- 
certainty in the estimate of the state vector is sufficiently small to allow 
the trajectory to be discreetly changed by adding the state vector to the 
position and velocity vectors on the integrated trajectory. This test is 
made to assure that large uncertainties in the estimate early in time will 
not produce the degenerate solution which is possible if noisy estimates of 
the state are allowed to be introduced directly into the system along with the 
desired results of the filtering and corrections for the non-linear nature 
of the problem. 


In order to produce this test, however, it is necessary to define a 
comparison function and specify a measure of error which is considered 
acceptable as a vhreshold for updating. This process, quite arbitrarily, 
has been performed as follows: 

3 
z Are y* avy] 

Fo: 2 [ ( 7) +> (=) 
and the weighting factor has been established so that the sensitivity of the 
semi-major axis to errors ir radius ard velocity for a circnlar orbit are cqual, 
Vis. 

4a,,2 2 4r22QAr 

a. ror br 


S298. 


ST) 65 1202-1 











fal a AV DV 
St 2(a-) He 2 


Thus, for this criteria, the two ratios in question should be weighted 
equallv, That is, X= 1. The threshold for this function has beer set at 


2 x 10” : 


Upon entry into UPSTAT, the covariance matrix for the errors in the 
state vector is examined; and the elements along th? principle diagonal (tr- 
trace) are extracted and utilized as Arn (or AV; ) in the construction °t 
the comparison function. The test is then made to determine whether updat i.,, 
will be accomplished; if not, the position and velocity vectors at this time 
are stored for future --~putation of the transition matrix, and wulil 
is returned to the calling routine. However, if updating is practical, it 
is performed, the conic reference trajectory is rectified, the state vector 
is zeroed, and the "classic" elements of the osculating conic trajectory 
are computed for reference (these variables are not utilized in the program 
due to problems of indeterminancies as e and w approach zero). Operation 
from this point is identical to the case where the test was failed. 


Computational Logic: 







UPSTAT 
write time, state vector, and 
estimation errors 


determine accuracy of estimate 





2320. 
STD 6& 1202~1 





Pied 








F>2x 10° 


ores 





write: 
there was no 
rectification 

at this time 








write: the trajectory has 
been rectified, 


CALL ELEMEN 





RETURN 
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2.4.1.5.1 Subroutine ELEMEN 


Purpose: to compute the classical conic elements (for the purpose 
of user convenience) at those times when estimation of 
the state vector has been performed (providing that 
tolerances established in UPSTAT for the estimation 
errors are satisfied). 


Deck Name: ELEM 


Calling Sequence: | CALL ELEMEN (R, V, GM, REQ, A, E, P, THETA, OMEGA, OMEGAW, 
OINC, OMDOT, WDOT) 


Input/Output : 


FORTRAN Math Common/ 
1/0 Name Name | Dimension | Argument Definition 
I R r 3 ARG 


radius vector (cartesian) Km | ' 





I Vv ¥ 3 ARG velocity vector (cartesian) 
I GM u 1 gravitational constant sii ‘ 
for central body (km3/sec2) 
I RQ Re iL equatorial radius (km) ; 
0 A a 1 semimajor axis (km) 
0 E e 2 eccentricity 
0 P 1 semilatus rectum (km) 
0 THETA r=) iL true anomaly of epoch (rad) 
0 OMEGA mee 1 right ascension of ascend~ 
ing node (rad) 
0 OMEGAW w L argument of periapse (rad) 
0 OINC a ee orbital inclination (rad) 
0 | OMDOT A 1 secular change in. due 
to earths oblateness (rad/rev) 
0 WDOT AW L secular change inw due 


to earths oblateness (rad/rev) 


~346- 


SID 65 1203-1 








Subroutines required: CROSS (cross product) 


: Functions required: AMAG (vector magnitude) 
DOT (dot product) 
ATAN (are tangent) 
SQRT {square root) 
SIN 
Approximate Deck 
Length: 64,0 (octal) 
Formulation: 


The paragraphs which follow present a summary of the equations which are 
mechanized in the routine and other information required to resolve ambiguities 
in quadrant, etc. This approach to the forwvlation assumes a degree of fam- 
iliarity with the equations of conic motion in terms of the "classic" 
elements; thus, should questions arise, a discussion similar to that of Ref. 

1 should be reviewed, 


The first step to be performed is the definition of the plane of motion. 
This set of computations will be accomplished through consideration of the 
angular momentum vector (per unit mass) and will assume that the position 
. and velocity vectors exist in a cartesian format in some selected frame of 
reference (the main program operates in the true equator of date frame). 


é hb=Fxy¥ 
H=h/{hl 
« The orientation of Hi can now be obtained in terms of the orbital 


inclination and the right ascension of the ascending node by two simple rota- 
tional transformations 
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Now equating the two expressions for # yields the desired information: 


COOwW = 4, 0242 </90 
ston Lf. = Mi fain 
sim € 


Coo S2 * Ua) 


The next step is the computation of the energy of the conic (or its 
equivalent the semimajor axis) and the eccentricity. This computation in- 
volves both the energy and angular momentum equations 


ee | VE ae a<O Hyperbolic motion 
SQ KR 

p={4] fis a?0 Elliptic motion 
e 


\! 
ps 


The final step is the computation of the argument of periapse ( W ) 
and the position in the orbit at the given time ( 9 = true anomaly of epoch) 
(see sketch) 





This evaluation is accomplished as follows: 


Y= coo"! E 2 i 


{r] 





Y= angle from ascending 
node =9@ *#® 
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The routine concludes with the computation of the secular perturbations 
in the right ascension of the ascending node (-2 ) and in the argument of 


periapse { w ) resulting from the first coefficient of the potential function 
( 


Ad - 37 J ( Be)" Cto x rad/rev 


" 


AW 37 Fi (Re)” (2-% anc) rad/rev 
P 


where: Re = the equatorial radius of the oblate spheroid 


References 


1) Townsend, G. E., "Orbital Mechanics" in the Orbital Flight Handbook, 
Volume 1, Part 1, Chapter 3 NASA - S®-33, (1963). 


Computational Logic: 
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W= p-e 
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2.4.1.6 Subroutine DATAPE 


Purpose: Reads a specially generated magnetic tape containing 
smoothed and ordered coordinate data, 


Deck Name: DAPE 
Calling Sequence: CALL DATAPE (KOUNT) 
Input/Output: 


FORTRAN anes Conmon/ 
1/0 | Name Dimension a) Description 


O 8 86©| TW The whole and fractional 
part of the number of days 
relative to 1950.0 at which 
the next observation will 
occur (1950.0 = J.D. 
24,33282.423) 






















0 ISTN The number for the observ~ 


ing station 
| O ITYPE The type of data: 
1, Range 

2, Range~Rate 

3, Azimuth, Elevation 

| 4, Range, Range-Rate 

5, Range, Azimuth, Elevation 
6, Range~Rate, Azimuth, 
Elevation 
0 ODATA The components of the 
observation vector 


Control indicator. 


I/O | KOUNT 









i, Do not return next 
point (INPUT) 

2, Return next point 
(INPUT) 

3, No more data points 

on tape (OUTPUT) 
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Subroutines Required: None 

Functions Required: None 

Approximate Storage Required: 650 (octal) 

Restrictions: Requires that the input data tape to be 
generated by FS4~305A be mounted on logical 


tape drive unit 9. 


Nomenclature: 


Dimensions Description 


A Buffer array containing a 
single logical tape record, 
IFIRST First pass indicator 


IGP Logical record counter, 
IGP=1, NGPS 


NGPS | Number of logical data 


records on tape. 


NPERGP Number of points per logical 
record, excluding final record. 


NPREM Number of points per final 
record. 


NPGINT Buffer array point index, 
NPZINT=1, (NPERGP or NPREM) 


XJDREF Pregram reference Julian 
date. 





Method: 

The first CALL statement to this subroutine mechanizes the input data 
tape prepared by the preliminary processor program. An information record 
containing J.D.2433281.5 (unused in this program), the number of logical re- 
cords, and the number of data points* per record is read into the program. The 
first logical data record is then read into the buffer array A and, after 
extracting the first data point from the buffer, control is returned to the 
eslling routine. Subsequent CALL statements extract single points sequentially 
from the buffer. After the final point within the buffer has been extracted, 
the next logical record is read in and the procedure is repeated until all 
data has been read from the tape at which time the count index (KOUNT) is set 
equal to 3 and control is returned to the norm program. 


#A data point is defined as the ordered set TW, TF, ISTN, ITYPE, (DATA. 
¥HID (2h) 33282.423 
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2.5 Matrix and Vector Algebra Group 


The routines of this group are general purpose algebraic routines 
for vectors and matrices, and there are no inherent limitations in the 
formulations or codings which restrict application (though there are 
fixed dimensions in some cases which should be altered to provide expanded 
coverage). Accuracy is assured in those cases where single precision 
arithmetic may produce loss of significances by employing double precision. 
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2.5.1 Subroutine MATMPY (Matrix Multiplication) 


Purpose: MATMPY is designed tc multiply any two conformable single 
precision matrices (with less than 70 elements each) to 
obtain the single precision product. All operations 
interior to the routine are performed in dorble precision 
to control roundoff and ‘oss of significance. 


Deck Name: MX PY 
Calling Sequence: CALL MATMPY (C, I, K, D, K, d, CD) 


Input/Output: 


FORTRAN Math Common/ 
Name Name Dimension Argument Definition 


(rows) Array of numbers to be 
(columns ) used in the premulti- 
plication of D by C. 


(rows) Array of numbers to be 
(columns) premultiplied by the 
matrix C. 

(rows) Product array. 
(columns) 

Subroutines Required: Nore 

Functions Required: None 

Approximate Deck 


Length: 700 (decimal) 
1300 (octal) 
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2.5.2 


Subroutine MTINV (Matrix Inverse) 


Purpose: MTINV computes the inverse of a nonsingular square array 


(of up to 36 elements) using the theorem which states that 
if a sequence of row operations will reduce a matrix to 
the identity matrix, then the same series of operations 
performed on the identity matrix will produce the inverse. 
(Error control is maintained internal to the routine with 
double precision arithmetic though input and output are 
single precision). 


Deck Name: INV 


Calling Sequence: CALL MTINV (B, ES, N) 


Input/Output: 


FORTRAN | Math 
I/O Name Name Dimension 


I B B N XN Arg 
0 ES Bl] NXN Arg 
i N N i 







Definition 








The N X N array of numbers 
to be inverted. (single 
precision) 


The N X N inverse of B. 
(single precision) 


The dimension of B 


Subroutines Required: CHOOSE (Check for singular B) 


Functions Required: None 


Approximate Deck 


Length: 


500 (decimal) 
7L0 (octal) 
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2.5.2.1 Subroutine CHOOSE 


Purpose: CHOOSE is utilized in conjunction with MTINV to determine 
if the matrix (of order less than 6) is sufficiently non- 
singular to allow the inverse to be constructed without 
excessive numerical difficulty. 

Deck Name: CHSE 


Calling Sequence: CALL CHOOSE (A, E, M, N) 


Math Common/ 
Name Dimension Argument Definition 


Input/Output : 





FORTRAN 
Name 





1/0 





A Arg The array of numbers which 
is being reduced to the 
identity matrix. 

E Arg The inverse being con- 


structed from A. 


- Arg A row counter for the 
operations being performed 
on A. 

~ Arg Order of the square array 
A. 


Subroutine nequired: None 
Functions Required: None 
Approximate Deck 


Length: 200 (decimal) 
270 (octal) 
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2.3.5 Subroutine ADDMAT 


Purpose: ADDMAT is designed to add two arbitrary vectors or 
matrices of the same order 


Deck Name; ADDM 
Calling Sequence: CALL ADDIMAT (A, I, J, B, C) 


Input/Output: 


Common/ 
Dimension Argument Definition 


(rows ) Two arrays of 
(columns ) numbers which are to 
be added 


(rows ) Matrix containing 
(columns ) the sum of A and B 





Subroutines Required: None 
Functions Required: None 


Approximate Deck 
Length: 80 (decimal) 
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2.5.4 Subroutine SUBMAT 


Purpose: SUBMAT subtracts the I by J array of numbers B from 
another array (A) of the same order 


Deck Name: SUBM 
Calling Sequence: CALL SUBMAT (A, I, J, B, C) 


Input/Output: 


Common/ 
Dimension Argument Definition 


I (rows ) A is an arbitrary 

J (columns ) array to be reduced 
by a second array of 
the same dimension 


(B) 


I (rows) CZA-B 
J (colwnns ) 





Subroutine Required: None 
Functions Required: None 


Approximate Deck 
Length: 80 (decimal) 
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2.5.5 Subroutine TRANSP (Matrix Transportation) 


Purpose: TRANSP constructs the transpose of an arbitrary array of 
numbers 
Deck Name: TRSP 


Calling Sequence: CALL TRANSP (A, N, M, B) 


Input/Output: 


FORTRAN | Math Common/ 
Name Name Dimension Argument 


A A N (rows) 
M (columns) 
M (rows) 
N (columns) 


Subroutines Required: None 






oR 


Definition 






Array to be transposed, 









The array containing 
the transpose. 


Functions Required: None 
Approximate Deck 


Length: 80 (decimal) 
140 (octal) 


~386- 


SID 65-1203-1 









Tot 


J9Vd 


OLOONVAL 
O9OCONVUL 
OSOQINVSL 
OVUONVGL 
O£ JONVUL 
OC OONVUL 
GTUUNVYL 


S8/ec/TT 


(S)NJI 


SID 65-1203-1 


ON3 

NUNLAY 

Vi= C14) 
wWfil=f T od 
N‘T=1I 1 0G 
(N‘W)O4 (NANI VY NOISNIWIG 


(rc *t) 


(Q°w*nf VIGSNVUL JNILNANNANS 


ANAWSLVLS JIINOS 


Nao = dSul 
SOESd 


2387. 


NOTIV3074 NsAI 


AdAL NOYLVION 


9LT 


39Vd 


’ 


SID 65-1203-1 


"LETOD SI Iwid0 NI HION37 9930 
LEOIO vl I 
N33 NO1LV301 NII N33 NOTLVI01 N3I N33 
JINIGNOdS3SYYOD NST NSS 
% NOLLIBS IOISAS 
G371V2 SINILNOWENS 
€ NUTLIAS dSNVUL 
SLNIODd AYLN 
I T0090 1 
TOBWAS IdAL NOTLVIOT TOGWAS 3dAd NOLLVION IDEWAS 
SAJIGVIUVA WVd90Ud G3INOISNSWIGNN 
dSNVUL 3NILNOYENS 
dVW 39VUDLS dSUl ce 
S8/EZ/TI SOES3 “ 


2.5.6 Subroutine CROSS 


Purpose: 


Deck Name: 
Calling Sequence 


Input/Output: 


Subroutines Required: 





Functicns Required: 


Approximate Deck 
Length: 


CROSS constructs the vector product of two arbitrary 
vectors A and B (in the following order C =A X B) 


CROSS 


CALL CROSS (A, B, C) 


Dimension Definition 


The two vectors to 
be multiplied in the 
sense A XB 


The vector product 


ereareerear 


None 


None 


70 
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: 2.5.7 Subroutine DOT 
an 


Nona 
fad Purpose: DOT computes the scalar product of two 3-vectors 


Deck Name: Not 


Calling Sequence: DOT (A, B) 


Math Common/ 
Name Dimension Argument, Definition 
A 3 


Input/Output: 


The two 3--vectors to be 
utilized in constructing 
the scalar product. 


The scalar product. 





Subroutines Required: None 
Functions Required: None 


Approximate Deck 
Length: 50 (decimal) 
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2.5.8 Subroutine AMAG 


Purpose: AMAG constructs the scalar length of a vector 
Deck Name: AMAG 
Calling Sequence: AM): (A) 


Input/Output: 


Definition 


The 3~vector (in cartesian 


coordinates) for which the 
length is desired. 


Vector magnitude. 


Subroutines Required: None 


Functions Required: DOT (Scalar product) 
Approximate Deck 30 (decimal) 
Length: 
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2.6 General Purpose Math Group 


In addition to the functions and routines required for algebraic 
manipulation of vectors and matrices (as presented in the preceding 
discussions), there are several mathematical functions which are utilized 
in conjunction with the conic motion formulation (and the related state 
transition formulation) and in conjunction with the definition of azimuth 
of the vehicle relative to the tracking station which must be discussed. 
These routines, due to their general nature and the fact that they are 
employed at several points in the program, do not logically fall in one of 
the other main groups. They have, thus, been presented in a separate section 


on the following pages (no formulation or computational logic due to their 
simple form). 
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2.6.1 Subroutine SINH 


Purpose: SINH computes the hyperbolic sine of a single precisio. 
argument expressed in radians. 


Deck Name: SINH 
Calling Sequence: SINH (X) 
Input/Output : 
FORTRAN 
Name Definition 


Single precision argu- 


ment (expressed in radians) 
for which SNH is to ve 
evaluated. 


Hyperbolic sine. 





Subroutines Required: None 


Functions Required: None 


Approximate Deck 
Length: 100 (decimal) 
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4 2.6.2 Subroutine COSH 


ae Purpose: 


Deck Name: 

Calling Sequence: 

Input/Output: 
FORTRAN 

1/O Name 

Z X 

G CUSH 


COSH computes the hyperbolic cosine of a single precision 
argument expressed in radians 


COSH 

COSH (X) 

Math Common/ 

Name Dimension Argument Definition 

X 1 Arg Single precision argument 
(expressed in radians) 
for which COSH is to be 
evaluated. 

cosh 2: - Hyperbolic Cosine. 


Subroutines Required: None 


Functions Required: None 


Approximate Deck 
Length: 





100 (decimal) 
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2.6.2 Subroutine COSH 


Purpose: COSH computes the hyperbolic cosine of a single precision 
argument expressed in radians 


Deck Name: COSH 


Calling Sequence: COSH (X) 


Math Common/ 
Name Dimension Argument 


Arg 


Input/Output : 


FORTRAN 
Name 







Definition 











Single precision argument 
(expressed in radians) 
for which COSH is to be 
evaluated. 


Hyperbolic Cosine. 


Subroutines Required: None 


Functions Required: None 


Approximate Deck 
Length: 100 (decimal) 
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2.6.3 Subroutine ARKTNS 


Purpose: ARKINS computes the single precision arc tangent (when 
given the sine and cosine of the angle) and assigns the 
angle to the proper quadrant (-72e<7 or 0z.@ 2277 ) 

Deck Name: ATAM 

Calling Sequence: ARKTNS (N, X, Y) 


Input/Output: 


FORTRAN Math Common/ 
Name Name | Dimension | Argument 
N - 1 








Definition 





Fixed point number to 
identify the range of the 
are tangent 

N= 180 -180*9° « 180 
N = 360 0£4e 2 360 











Cosine of angle 
Sine of angle 


Arc tangent (Y/X) 


Subroutine Required: None 


Functions Required: ATAN (arc tangent function) 


Approximate Deck 
Length: 160 (decimal) 
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2.6.4 Subroutine DERAQ 


Pa 
So ee 


< Purpose: DERAQ is intended to represent the Kronecker delta 
ae (i.e., bi ) and as such has the value of 1.0 or 0.0 
ne depending on the arguments I and J 


Deck Name: DERQ 


Calling Sequence: DERAQ (I, J) 


Input/Output : 







FORTRAN | Math Common/ 
Name Name Dizension Argument Definition 


Lied des, ake Two fixed point variables 
defining the value of. 
DERAQ Kronecker delta, 


Subroutine Required: None 









Functions Required: None 


Approximate Deck 
Length: 30 (decimal) 
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3.0 Accuracy Tests 


The program logic has been checked numerically against three precision 
programs, 


1) AP1LO - A general purpose trajectory program prepared for the 
Apollo program and employing an Incke formulation and an 
Adams integration logic. 


2) SPACE - A JPL single precision cowell trajectory program employing 
a kth order Runge-Kutta integration logic. 


3) ITEM - A Goddard trajectory program employing an Encke formulation. 


The results of these tests indicete that agreement can always be obtained 
through the sixth digit (frequently through the seventh) for integrations 
covering several orbits of the earth providing that the constants of the 
programs are comnatible. This level of agreement is considered as a proof 
of the program logic. This conclusion is strengthened when it is observed 
that the basic logics of the three programs are somewhat different and that 
in no case could the agreement be expected through the eighth digit due to 
the restricted word length employed in single precision on the IBM 7094. 
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4.0 SAMPLE PROBLEM 


4,1 SAMPLE PROBLEM DISCUSSION 


Punched paper tapes with tracking data recorded by the Floyd Terminal for 
passes 6069 and 6070 of ECHO II were provided for the purpose of demonstrating 
the operation and accuracy of the differential corrections program described 
on the previous pages. Thus, the discussions which follow will enumerate 

the procedures involved in construction and solution of this sample problem. 


4.2 ORBIT PARAMETERS 


Passes 6069 and 6070 of ECHO II occurred on 27 April 1965. Thus, published 
data for this satellite was searched to determine an estimate of the position 
and velocity at some epoch prior to this date. This procedure was 

required to assure that the trajectory could be generated accurately so that 
the differential corrections process can converge to the desired trajectory. 


The Goddard Space Flight Center on 30 April 1965 issued the following 


information pertaining to 1964 O4-A (ECHO II) in TWX NR 054/301-474-4911 (the 
true equator of data frame of reference). 


Epoch OMS UT, 20 April 1965 


semima jor axis 7528.31 Km 
eccentricity 0.02447 
inclination 81 .4.50° 
Mean Anomaly 113 .092° 
argument of periapse 34.156° 
right ascension of 31.202° 


ascending mode 


Accordingly, this information was resolved into the corresponding position 
and velocity vectors. (Using subroutine POSVEL) the resvlt was: 


Epoch Oo% uT, 20 April 1965 


Radius Vector = -5915.9438 - 2918.1582  3783.0742 (km) 
Velocity Vector = -2.7144510 ~-2.7299969 » - 6.0748353 (km/sec) 


At this point, this information was updated to an epoch just prior to 6069 

by employing the general perturbations program described in SID 1203-3, 

(The elements themselves could have been updated employing only the secular 
changes inw and.%. However, this approach was discarded to avoid introducing 
the slight errors associated with such a process.) The results of this 

process are: 


Epoch 27.63865740, April 1965 


Radius Vector = 4952.3943 1406.9609 -5362.9226 Km 
Velocity Vector = 4.4573218 2.9062537 5.0928345 Km/sec 
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4.3 REDUCTION OF THE RAW DATA 
4.3.1 RAW DATA TAPE 


Input deta for the sample problem consists of two paper tapes recorded 

by the Floyd tracking facility containing the following observations for 

the ECHO II satellite passes 6069 and 6070; elevation and azimuth in degrees, 
doppler reading in kilo-cycles, a lock-on indicator, and universal time in 
seconds. Since the IBM 7094 computing system utilized does not have a 

paper tave input capability, a short IBM U,01 program was written to transfer 
the data directly to magnetic tape with the format shown in Figure 1 (one 
physical file per pass). The leading information record for each data file 
contains the following indicators. 


PASS PASS 
SYMBOL DESCRIPTION 6069 6070 
NSTN Code indicating the tracking a8 1 
station recording the data 
(Floyd, see input) 
NTYPE Code Indicating the type of 6 6 
observed data (range rate, 
azimuth, elevation) 
NNUM Total number of words per file. 3561, 3392 
(his was determined by tne 1,01 
program in transcribing the 
data from punched paper tape) 
XJDATA Julian date (zero hour UT) of (24 )38877.5 (24) 38877.5 


the first observation for each 
file. (27 April, 1965) 


A third file containing NSTIN equal to zero indicates to the program that 
all data has been read from the tape. 


Graphical representation of the raw data is presented in Figure 2 

through T (the broad line apparent in these figures is in reality a 

series of points or plotting symbols). Examination of Figure 6 will dis- 
close random irregularities of azimuth observations being recorded 360 
degrees out of phase, especially prevalent for angles approaching 360 degrees 
(e.g., +354° recorded as -6°). These points were adjusted in the preliminary 
processor inmediately prior to smoothing the data segment to assure that a 
consistent convention is employed in recording the observed and computed 
values of the data. 
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TRACKING STATION 

TYPE OF DATA 

NUMGER OF WORDS PER FILE 
JULIAN DATE OF FIRST DATA POINT 
TIME 


CGORDINATES DATA 
AS INDICATED POINT 
BY NTYPE 


NSTN20 INDICATES NO MORE 
DATA ON THIS TAPE. 


* FORTRAN & 

¢ BINARY MOOE 

* MULTI- PHYSICAL FILES 
‘LOGICALLY PACKED TIME 
AND COOR DIMATE DATA 


STD 65-1203-1 











This figure also serves to illustrate a data format which is net generally 
acceptable to the program. Although most of the data are observed to be 
sequenced chronologically, a few recording errors of the time are evident. 
At this time no attempt has been made to determine if such errors exist. 
Rather, it is considered to be the user's responsibility to assure that 
such errors have been eliminated. 


Fortunately for this problem these errors do not affect the solution due 

to the manner employed in smoothing the raw data over short intervais of 
time. However, this circumstance exists only because the erroneous time 
values did not occur at the midpoint of a 20 consecutive point strip of 

data (see subroutine FIT in the preprocessor SID 65~1202~2). Since such 
fortuitous arrangement of errors cannot be guaranteed, remedial steps should 
be taken in recording the data. 
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4.3.2 Variable Dimensions for Preprocessor 


The required values for the variable dimensiors were determined by the 
procedure outlined below. (See "Variable Dimensions" and "Storage Limita- 
tions" in the "Program Operation" section of SID 65-1203~2) 


AA(4, MAXAA): This array must be large enough to hold the smoothed roints 
corresponding te a single raw data file. The largest raw data file contains 
3564, words. The number of words per information record is 4, and there are 
4 words per data point (time, doppler reading, azimuth, elevation.) Thus, 
the number of raw data points per file is (3564-4)/4 = 890. Now the 
smoothing routine reduces 20 raw data points to a single smoothed point, 
therefore, since 890/20 45, MAXAA was dimensioned by 50. 


A(6, MAXA): The primary purpose of this array is to receive the smoothed 
data from the AA array. Since MAXAA was set to 50 and there are two files 

of data, MAXA was set to 2x 50 =100. If there had been storage limitations, 
MAXA could have been reduced (the smoothed points would then be temporarily 
stored on tape). However, this array must be large enough to hold the 
smoothed points corresponding to a single raw data file. 


STN (6, MAXSTN): The usual criterion for dimensioning this array is number 
of raw data points per file. The sample problem has MAXSTN = 1000. 


This procedure is outlined here since more demanding tasks will undoubtedly 
be presex.ced to the program during its period of use. The sample, in no 
way taxes the capability of the program. Before passing it is noted that 
these variables can be assigned large values to assure operation for most 
problems and that should these limits ever be violated a message will be 
printed as to the nature of the problem and steps required to correct it. 


4.3.3 Input Data Load Sheet for the Preliminary Processor 
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4.3.4 Preprocessor Output 


Primary output of the preproc«ssor is a magnetic tape read directly by the 
differential corrections program (a format shown in Figure 8), However, 
several print options exist (see SID 65~1203-2). These options include 


printing of the raw data 
printing of smoothed data 


Since the raw data were provided in graphical form on , zevious pages, the 
irst option was not mechanized. The second option was however employed. 


Accordingly, smoothed and sorted data corresponding to the curves rresented 
earlier are presented in the following table. 


a Fe) 
SID 65-1203-1 


en tt tne 


FIGURE 8 


FS4-305A 
OUTPUT TAPE FORMAT 
SMOOTHED AND ORDERED DATA 





XJDREF 
NGPS 
NPREM 


NSTN 







PROGRAM REFERENCE DATE 
INFORMATION 


TOTAL NUMBER OF RECORDS 
RECORD 


NUMBER OF POINTS PER RECORD 

CLUDING FINAL RECORD 
NUME . OF POINTS IN FINAL RECORD 
WHOLE DAYS FROM (950.0 
FRACTIONAL) JO. 2VIS2f4.7AS 
STATION 


TYPE OF DATA DATA POINT 


DATA RECORD CONTAINING: 
*NPERGP POINTS 
* SEVEN WORDS PER POINT 


AS INDICATED 


ks inoieateo | 
BY NIYPE 


* FORTRAN [Wy 
* BINARY MODE 


* SINGLE FILE (NO FILE MARKS) 


~h25- 


SID 65-1203-2 


00 JTELT606E2°O 
00 318¢2T9E9C"O 
OO 368292062°0 
00 J68008ETE*O 
00 3J1SSO96EL°O 
OG 3ALZ90LY9IE*O 
OO 3S802¥68E°O 
CO 3090S6ET?°O 
00 3€28ECSE4°O 
06 316896%S7°9 
00 3TT6LYILY°O 
00 3¥7099098%7°0 
00 389S¢TS64¥°O 
00 3099T0S64°0O 
00 36¢0%79726%7°0 
00 38€062S84°0O 
00 390886¢cL4°0 
00 3¢28S4%18S7°0 
00 3926€89t7°0 
00 39289¢4¥T"O 
0O 3S6T2L68E°O 
00 3S2062%9€°0 
00 3100689EE°O 
00 3svecc¥yITe’O 
00 369261S8C°0 
00 36T08TESZ°0 
QO 3676SStt2°0 
00 3TECOEsso7e*o 
00 3€926S498t°O 
00 3679C¥7T9T°O 
OO 396C6TSET°PO 
00 3ZLT909TI°O 
T0-392800S66°0 
T0-36€Sc00SL°0 
TO-30€S%796¢SS°0 
TO-3T6649S41E°9 
1O-31cascyo7 "°C 


Z 


0O ALETB8bT9°O TO 398000884°0 9 T 00 309S299¢L °0 
00 3261S591569°0 ~ “TO FO8SZ669F*0 FT 00 SEEHYYIEL*O 
00 34%029%7002°0 TO 3¥9EBESYH"O 9 T 00 36S2Tz29EL °O 
00 JSTO6O0LHL°O “10 3S6TTESZTH* OF TT 00 SETSESEL*O 
00 36S596€862°0 TO 39062200%7°0 9 Tt 00 38S6%7LSEL°O 


00 SHOECESSB"O I0 SOLLE90LE*O ~ 9 TT 00 JTEBSTSSEL*O ~~“ ¥O 30000S6SS°0 


0O 31220SlT6°0 


10 380ZESlee*O 


9 


Tt 


00 


3245998 2S¢2 °0 


00 391z7zesEE"O ~—" TO 3LE6cTHOUE"O CT 0 TE SSOSEL *O 


TO 3e€v0c090T°O 10 3TLTO86S2°0 9 T 00 AvOyYCBYEL °O 
TO JO0ESZOVTT°O TO JBecyeste"0 FT TT ~ O00 30ECESPEL*O 
10 30c68S¢cT°O TO 3S6€9CL97T°O 9 T 00 A£OT9EVELTPO 
IO’ JIOORITETTO ~ TO SLERHvOTT "O_O 00 66 CTHEL OO 
TO JOL9SOTYT°O 00 JS9LEVYVESD°O 9 T 0O 320868Eeel°O 
TO 3S4¥€690ST°O 00 SS6VIS9NSTTO FF ~~ T ~ 00 38¢999€E2 °0 
10 39LSLE09T°O 00 3980896479 °0- 9 T 00 310SEyvEEL °O 
10 3S26¢2669T°O ‘TO SLTLO660T*O- = —s 9 tT 00 SicEeOZzEeEL °O 
TO ALESSTELT°O TO 30SZ2ECLST°O- 9 T QO 300¢262EL*°0 
IO STILLS6L8T°0 ~~ TO SEOTIOSOSY*0="~— GT OOS YLOFL CEL *O 
TO J0CTS79ET°O TO 389668642 ° 0- 9 T 00 36680SceL°O 
TO 3928L6E02°0- TO 3LSL98%62°0- 9 T O00 3€lLl2cel °0 
10 3€0790TTZ°0 10 AEvELOTEE*O- 2 Tt 00 366S6470¢2E1°0 
TO 3SET6SLTC°O TO 3tT200L9ggc*0- ~— 9 : 00 32l¥Tstel*o 
TO 38L1¢SE722°0 TO Jils9sSB6E *0~- 9 T 00 386 c8sTttd °O 
10 3212768¢7°9 TO J€664L924°0- “Gg ~"T 00 3Tz2IsEeTtel*o 
10 392¢988EeEC °O TO 300T¥22S%*° 0- 4 t OO 32661TItL°O 
TO 3LE9¥VERES? PO TO AVETLYYLY*° O-— 9 T 00 30L8880EL°O 
TO 3209e4249C*°0 10 382020%64°0- 9 T 00 3969S90€L°0 
TO JOS8ET9~C°O TO 3€2980¢cTS*0- 9 T 00 369Sc70tL °0 
TO 39LBCS64HC°O TO 36%97¢L12S°0~ g T 00 sevyeTOEl °0 
TO 3499492S2°O TO 300%SS0%S °0- 9 T 00 389296621 °0 
TO 36699%SS¢°0 TO 36¢cSL61SS°0- 9 T OO 3cvtcl6eLl°O 
TO 362€S28Sc°O TO 3LG6668729S°0- 9 T 00 3E7e8e76el °O 
TO 3862199092°0 TO 30S6266TLS°0- ce] T 00 3679S262L°0 
TO 366606292°0 TO 3612472626 °0- 9 Tt 00 32¢26206¢22°0 
TO 329¢¢C6492°0 TO 34¥02998S°0- 9 t 00 396€62822°0 
TO 312288992 °0 TO STLL8ECES* 0O- 9 T 00 312¢9Sbel °O0 
TO 342822892°0 TO 3471¥€916S9°0- 9 T 00 3S60CE82L°O 
7 VES 
A x ddAlL NLS 


eze*cszecve SI 3LVG 329N3¥34354 


see VWLVG G3LdOS GNV GSHLOUWS xxx 


4O 30000S6SS°0 
~ ¥0~ 39000S5SS°0 
70 30000S6SS°0 
40° 30000S6SS°0 
40 39000S69S°0 
7Q 30000S6SS°0 
70 30000S6SS°0 
70 30000S6SS°0 
70 30300S6SS°"0 
¥O 30000S6SS°0 
90 30000S6SS°0 
¥O 30000S6SS°0 
70 30900S6SS°0 
70 30900056SS°0 
70 300609S55SS’0 
70 30000S55S°0 
70 30000S5sS5°0 
70 30000S5SS°0 
70 30000S6SS°0 
40 30000S5SS°0 
70 30000S6SS°0 
70 30000S6SS°0 
70 300V0S6SS°0 
70 30000S6SS°0 
70 30030S6SS°0 
70 39000S6SS°0 
70 300006655 °0 
70 30000S5SS°0 
40 30000S6SS°0 
%O 30000S5SS°0 
7Q 30000S6SS°0 
¥O 3090056S3°0 
70 3900056555 °0 
4O 3000055S%5°0 
70 30000S65S°0 
70 30000S5SS°0 
J 1IDHM 


(3iVG 35x WOXA SAVO) AWIL 


LE 
9€- 
Ge 
HE 
ce 


--9¢ - 


SFHNMTN OK DRMOHNMT MOF OBAO er 
BAA MMA MH NN NN NA AOL AN OM 
SID 65~-1203-1 


ANOS Or CEO 
ed 


ON 
INI 


426 





4 
1 
' 


we we ee 


OO sccv6Tt8S°O 
00 329968T29°O 
00 3L1T2OT99°0O 
00 4¢S799S69°0 
09 3Te¥v1S72l°O 
00 32989¢¢7L°0 
00 3S0SE8691°9O 
OO 301c00S%1°0 
00 3i2l¢2982el°O 
QO AVTT98E0L°O 
00 3826T18199°0 
00 3469L4829°0O 
00 s3¥ySEY9BB8S°O 
00 32€€66247S°O 
00 36L09L€67°0O 
00 3S0€92bS4°0 
00 3BIySesty’°Oo 
00 3aLv6elb6LE°O 
00 3”v8EervvEXE°’O 
QO 3SS99SEOEe’O 
00 398472249L2°0O 
00 Aalst Ilsy2°O 
0O A6E€TES9ITC°O 
00 30T99T6BT°O 
00 399TO099ST°O 
00 36T2z08ET°O 
00 3XST9S+FTT"O 
10-3600 4881T6°0O 
TO-3¢ 40 TH 189 °O 
TO-3696478EL°O 
TO-3AL8E8c028°0O 
00 38SL6E90T°O 
00 3ELlEexvB9cT°O 
OO 3TETC6LYT°O 
00 30286cELT°O 
CO 3cLTLTC6T°O 
00 3as98¥vESTC*O 


Z 


ecy°*ceszeceve S1 JLVO 39N3uN3538 


TO 
TO 
10 
TO 
TO 
TO 
10 
TO 
TO 
TO 
10 
TO 
1s¢) 
TO 
10 
tO 
TO 
TO 
TO 
TO 
10 
TO 
TO 
TO 
TO 
TO 
TO 
TO 
TO 
0G 
00 
00 
00 
00 
00 
00 
00 


3L0880S25°0 TO 3CHYETSSE°O 9 00 39SLOBOT8°O 70 309909S6SS°0 
3L660299S°0 .0O JS6EO¥STIEC*O 9 : 00 36¢9LSO0T8°O %O 39090S5SS°0 
d22tS09S59°0 10 398¢28192°0 9 ( 00 3SS77E0T8 °O 70 J0J00S5SS°0 
J3e6C0S677S°0 'O 266¥719S02°0 9 { 00 36cETIOI8°O 70 39000S5SS°0 
3LV7ELSTES°O .O A8’Be00sSyT°O 2 T 00 3SST88608 °0 70 3900066SS°0 
JOvsovlTS*O 00 300494861 °C ? T 0O 4820S9698 °0 70 30000S6SS°0 
321S1TL20S°0 90 JOO0cELSct°O ? T 0O 37S8T7608°O 7Q 309000S6SS°0 
JTTSOBL&EY°O 90 3STS¢@yT 22° O- 7 T 0O 322248 1608°O %0 39900S6SS°0 
SOUZ6ceL9°O LO J3ESSOSOET’°O- ? T 00 3009G6808 °0 99 39990%6S59°0 
38S9166S9°9O TO 34¥90SS06T°0O- ? T 00 3927¢18908 °0 70 39000S6SS°0 
JLISSeiyy7°O TO ATLLY¥~CBH?C°O- 3 T 0O 30067898 °O ’O 39900S6SS°0 
AIS9CSTLEV"O TO 3TL6BOTOE*O- 2 T 0O 352192808 °0 ¥O 30000S6SS°0 
A8O9¥ER~H°O 10 39€2619¥VE° 0- 9 T 00 366620808 °0 760 30020S6S5°0 
J0Zvo9619°O TO 3¥5L¥S06E °O- 9 I 00 3728612108 °O 70 39000S6SS°0 
3Z2L969¢14¥°O TO 371T0C8924"°0- 9 T 00 38699S108°0O 70 39900S6S5°0 
39712 £990%°0 TO JAIiLyEORSY° O- 9 t 00 34%¢2SECL08 °O 70 320000S6SS°0 
3BeS¥vETO¥’O 10 34u9L¢284984°0- 9 T 00 ALECOTLOB°O 4O 39000S6SS°0 
3AZ4HL996E°O TO 394S47810%9°0- 9 T 00 3ECCLYEIOB°O 70 33900S5SS°0 
3€9S0926E°9 TO 39T89G622S°0- 9 T 00 4960479908 °0 70 34000055SS°0 
32426688C °O TO 3S211T¢274S°*0- 9 T 00 369609908 °0 70 30000S6SS°0 
3Z699LS8E°O TO JI2ve2sess*o- 2 if 00 396111908 °0 70 30000S6SS°0 
JTILE26CBCPO TO 3¥8€S¢TLS* O- 9 T 00 3699%76S08 °0 7Q 3990056SS°9 
3TcOce O8E°O TO 3cT8edcs8s°o- 9 i 00 37671TLS08 °O 70 3O000S6SS°0 
JeCCOOBLE °O TO 3”v¥202T6S°0- 9 I 00 389£847S08°0 70 30000S6SS°0 
39G6S8S2€°O TO JILLEBBES°O- 9 T 00 3€61S72S08°0 7G 30000S6SS°0 
3162 T6ele°O TO 3¢¢c€¢c4S09° 0- 9 T 00 349020S08 °0O 70 30900S6SS°0 
J698¥T CLE PO TO 30S9TETIO°O— 9 T 00 326882408 °0O 70 30000S6SS°0 
38640S 01¢€ °0 10 3LT6909T9°0- Le | 00 3991554908 °0 ¥O 30900S6S5°* 
SVECOLEIE °O TO 3€0€040C9°0- 2 T 00 3e0COEVO8 °O 70 30000S6SS°0 
JTLYYCET YH’ O 10 3900666L4°0O _9 ‘TT 00 3208S4%8€2°O 70 39000S6SS5°0 
4S61869¢%7°O LO 3LLS7BEIS*O 9 T 00 368S6c8el °O 70 393990S6SS% °0 
ASGLE66L94°0 TO 362LTELSS*O a T 00 329¥908EL°O 70 30000S6SS°0 
J€860L014°0O 10 39€9508%9S°0 9 T 00 3288cec8le2 °0 70 30000S6SS°0 
389S8%7S69°0 TO J8Sbe6EES*O ~ i T 00 3291092¢2°0O 70 390900S6SS°0 
ASLESECCS°O TQ 38Lt0L62S°0 9 T 00 sscO/Elel °O 70 30000S6¢SS°0 
B6L9641SS*O0 “TO JZHSBOLTS*O FF T OO ST98ETLEL*O™ — 0 399005555°0 
3TS8S9E8BS°O TO 3c¢c6Lyeos*d I 00 3vELO69EL°O 70 30000S5SS°0 
3VUs 3 IOHM 
Fe eae ae “Sd AL Nis —Tsiva 33a wWOUS SAVG) AWIi 
"ae Vivo GSLgoS GNV GSHIDOWS 44% 


Lh 
el 
ed 
12 
OL 
69 
89 
19 
39 
S9 
%9 
£9 
29 
9 
09 
6S 
BS 


-4.27- 


SID 65-1203-1 


3S 
SS 
%S 
€S 
cS 
TS 
OS 
64 
87 


97 
SY 
vad 
€¥ 
cx 
TY? 
Ov 
6€ 
8 


ON 
INI 7 


me 


Ri Sa ta i el anh tn or eee hell 


oe Ten > we. § 


SID 65~1203-1 


a rareibeieammneeenesatneaensattiilt ema tumemeeamienmneeneemetaes cate 


T9000 J3SGWNN dwng 


00 360€c80ST°O TO 3té2c0elce9°0 TO 3€eclsSo0se’o 9 T 00 369TS8ET8°O 7O 39990S5SS°0 68 
00 374297¢TST°O TO 36629T129°0 TO 36SLTOTL¥°O 9 T OO 3ScOvsEeTs°O 70 30000S6SS°0 8g 
00 3¢9S88C9T°O TO 3¥SETO9C9°O TO 3L26222115°0 9 T 00 JC6CTLETS°O 70 329030S6SS°0 ie 
00 398cEevsst°o TO 3S9620429°0 TO 388969€LS°0 ee | 00 J6TIS¥ETts°O 70 39090S6S5°0 98 
OO 3182S2St2e°o TO S30¢cT28T29°0 TO 3S€tE€6L9S°0 9 T 00 3266%2ET8°O ¥O 3000056Ss°0 S8 
O00 3STZTIT¥?*O TO 3SS9T¥6T9°O TO 3L44846SS°O 9 JT 00 38T8toEets°o 70 30390S6SS°0 98 
00 320619692°0 TO 3€89Sl9T9°0 TO 3811096%S°0 9 T 00 316981218 °0 7O 39090S6S5°0 £8 
00 360S5S966¢°"0 TO 38v7ETBCT9I"O TO 3LS2LlLES°O 9 ““T™ “00° 321959218 °0 7O 39000S5SS°0 28 
00 39€6800E£°O TO JECSLIOII"°O TO 3€2S66S92S°0 9 T 00 JT6ECEZIS°O 70 30030S5SS°0 Ts 
00 390646S9E°0 TO 3Te9EECCI"O TO 389S4260S°0 9 T 00 397T260c2T8°0 7O 309900S5SS°0 08 
OO 311S6200%°0 TO 3S6¢2SOIT9°O TO 38C99868¥°O 9 T 00 306098TT8°O 70 30029S6SS°0 6L 
00 3¥O0EBSLEY°O TO 3€S¢26S66S"0 TO 399TYCB9IY7O 9 t 0O 3STtS6c9TTts°0 ¥O 30000S6SS°0 BL 
00 328259919°0 TO 3479899¢6S°0 10 36¥76SS¢%9%7°0 9 T 00 36826ETT8°0 70 30000S5SS°0 dl 
OO 3888Cc89TS’O TO S6G02S98S°O° TO S8OLCESTH*O OG CT 00 32999TTT8°0 70 30000S6SS°0 9L 
00 3c0€S909S°0 TO 3TT?@vLe8S°0 TO 36SO09LE6E°0 9 T 00 aTe6cOTtis’o ¥O 3099055S3°0 SL 
Vas a IOHM ON 
Z A ° Xx ~AdAL NiS 7” (31VG sau Weds SAVG) JWIL 3NI) 


ecev°’cBceese SIT 3LiVd JINJdsd3a eax VIVO GALN0S ONV GSHIOONS «4% 


108 | 


PERT 


LZ —— 


SSamtnmtiereds MEI alter ye 
nda 


L.4 Sample Differential Correction Program Input 


The input for the sample problem wil] be discussed in relation to the 
significance of the data and its location in the fundamenta! data array. 
All quantities which can be readily changed wil] be discussed though some, 
as will be noted, will. not be changed from their preprogrammed values. 


Reference to the man of the fundamental array (Appendix 1 ) will varify 
the DATA locations given in the following pages. 
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h.4.1 Physical Constants for Math Model 


Equatorial Radius 
Polar Radius 
Coefficient J 

H 

D 
GM for the Earth 


Spin Rate for the Earth 


GM for the Moon 
GM for the Sun 
The astronomical Unit 


DATA (1) 
DATA (2) 
DATA (3) 
DATA (1) 
DATA (5) 
pata (6) 
para (7) 
DATA (8) 
DATA (9) 
DATA (10) 


As Programmed 
As Programmed 
As Programmed 
As Programmed 
As Programmed 
As Progranmed 
As Programmed 
As Programmed 
As Programmed 
As Programmed 


4.4.2 Satellite Data for ECHO II 
Mass (assumed) DATA (16) 1.0 
Area ecege DATA (17) 100.0 
Drag Coefficient (") DATA (18) 2.0 
Reflectivity (") DATA (19) 1.0 
Position Vector X DATA (20) 1952.3943 
¥ DATA (21) 1406 .9609 
Z DATA (22) ~5362.9226 
Velocity Vector X DATA (23) 4.4,5732).8 
Y DATA (24) 2.9062537 
Z Pata (25) 5 .0928345 
Time Relative to 1950.0 
OM April 27 1965 = 2h38877.5 
Epoch = ___.63€ E5740 
J.D. of Epoch 2438878 .13E 65740 
Reference Epoch  2/,33282.423 
299571565740 days 
TW DATA (26) 5595. days 
TF DATA (27) -71565740 days 
14.4.3 Control Options 
WINDEX DATA (28 1. 
CHECK DATA (29 l. 
GONO DATA (30 l. 
CODUMP DATA (31 l. 


4.4.4 Station Identification Card (SIc) 


Station one (Floyd) will be utilized 
Stations two through ten will be deleted 
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4.4.5 Station Data™ 


Latitude 


(rad) DATA (36) 75393225 
Longitude (rad) DATA (37) 1,.96824,718 
Altitude (Km) DATA (38) A797 
Station Name DATA (39) As Programmed 
Horizon Correction DATA (77) As Programmed 
Variance in Latitude DATA (87) As Programmed 
Longitude DAPA (88) As Programmed 
Altitude DATA (89) As Programmed 
Range DATA (176) As Progranmed 
Range~rate DATA (177) Ol HH 
Azimuth DATA (178) As Progrewed 
Elevation DATA (179) As Progroumed 


4.4.6 COVARIANCE MATRIX FOR ERRORS IN r AND VW 


The data for the elements imply that estimates of the errors in these para~ 
meters are approximately: 


36,, = -1 Km 
36%¢ = -0001 
.Ol rad 
30,m “ 5 rad 
36; 0.01 rad 
Aw - oe 
5163 
3a = a rad 
or eA 
. -9 
O,, ix lo 
gy. = (35x 10° 
Ac use 
rape e .35 x 1078 
x. an an~ 
Taw 2 0D X LV 


* See subroutine INPUT and Appendix 1. 

#4 From conversations with RADC personnel it was learned that the 
data from which R was computed were not as precise as had been 
earlier. No firm estimate of variance was available, however. 


Doppler 
assumed 
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Now, since this set of errors is relatable directly to the set of errors 
dy and dy in rotating coordinates by the matrix equation of Table 1 , 
and since the rotating and inertial coordinates are relatable by a linear 
transformation, the covariance matrix for d? and d¥ can be evaluated as 


{ax} al ae } 
{ ax’} 7 { axt = TA {ae} 


E(ax’ax’’)= TA E( AEE’) a'r’ 


tt 


u 


where AX = { ; ry rotating coordinates 
V 
dr ; 
Ls a cartesian coordinates 
A * matrix of Table 1 
aa 
AE = { Ae } 
si . Transformation from rotating to inertial coordinates 
Q3 
But ECAx’A x’ “) is not the matrix required because it is valid for 


the initial epoch rather than the updated epoch. Thus, one final transformation 
is required 


AX, * @(tgst,) AX, 
and 
é oT 7. 
Eca yy 03°) = Pl tt) E(ax/Ax, ) p (eyt, ) 
where y(t,; ¢,) is the state transition matrix. For the present purposes 


(tz 5%,) can be assumed to be adequately approximated using a conic reference 
trajectory (see subroutine TRANS). 
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| Adopting the notation T.(o4) etc. to mean a positive rotation about 
4 the X axis through the angle “«, the transformation of position errors 
becomes 


Boe T6-@) 7, Ge) Z[-Ce@+’%)] df 


2 R(2.,4,;¢@)d?F 


dv‘ = B-n) Hc~-2) GZ [Co + wy] & ("90 49) dv 


iM 


Reisen es) Ow 


where 


— 
ei 


= -G0 -(9+tw-—y) 


thus, T becomes 


Rle,4£,¢9) O 


0 R(e@,L, =) 
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These data must be re.d into the differential corrections program column - 
wise. (The digital computer stores elements of Ay 5 in a single subscripted 
array Qiny 214 93) 5 Qny) Qi24 Oa2y O32 -:- ane 3 etc.), thus 


Py = DATA (355) 
Po; = DATA (356) 
P37 = DATA (357) 
Py = DATA (358) 
Pe, = DATA (359) 
PZ) = DATA (360) 
P12 = DATA (361) 


4.4.7 Input Data Load Sheets for the Sample Problem 


As was noted in the discussion of subroutine INPUT, most of the data required 
for program operation can be prerecorded (in fact data for a check network 

is presently provided within the program). As was also noted in the discussion 
of INPUT and REED, these data will not be affected by inputting the desired 
date. unless new data are assigned to the locations presently occupied by 

these numbers in the common array (see Appendix 1). Thus, only the following 
data are required (the load sheets are in the proper format and should be 
copied) - 
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Five notes are worthy of attention before this discussion ends. 


1) 
2) 
3) 


h) 


5) 


The SIC must be the first physical card in the data deck. 
The 999 control. card must follow all numerical data. 


The 9999 control card must follow the station name data even if 
there are no changes (as in the case of the sample). 


The order of the cards following the SIC and preceding the 999 card 
can be arbitrary since a location is provided with each group of data. 


With the exception of the locations and data provided on the SIC, all 
numbers are floating point. These floating point numbers can be 
expressed with a decimal point (as on card 2) or without (as on 

cards 6-13). However, in the latter case, a decimal is assumed 

to be located between the first two locations of the field (i.e. 
before the first digit) and the last 3 locations of the field must 
contain the exponent of 10 which properly defines the number in 
question. 
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4.4.8 Program Output 


The differential corrections program was utilized in conjunction with 
the data of the preceding pages to compute the orbit of ECHO II. The resulis 
of this effort will be presented in their entirety to facilitate checkout on 
other systems and to demonstrate program operation. 


However, before presenting these data, it is felt advisable to discuss 
several problems which were encountered. First, the trajectory,as computed 
from the data provided ,checked the Goddard estimate of this time with negli- 
gible error and the observed values of azimuth and elevation to within 
one-half of a degree after seven days of prediction. This fact implied again 
a good level of accuracy. But, the computed values of range-rate failed to 
show this level of agreement. This fact was not considered entirely un- 
reasonable at the time. 


Second, reference to the raw data disclosed that the range-~rate as 
recorded exhibited a discontinuity of approximately 2.6 Ke at the time the 
doppler changed signs. This jump corresponds to approximately 0.172 Km/sec, 
and thus could explain the range-rate disagreement which was experienced if 
it could be ascertained how this correction should be applied; (i.e., is the 
positive segment of doppler correct, the negative correct, or are both in 
error?). Regardless, however, the data are in error and will affect the results 
in two distinct ways: 


1) Since the data are biased, the filter process as coded,is not 
correct (biases are assumed to have been subtracted from the 
observed data). 


2) Since the data are to be processed with the biases included, 
an oscillatory nature of the solution is to be expected (under 
the assumption that neither branch of the doppler curve is 
correct). 


Third, since the range-rate data was assumed to be in error, a call was 
placed to RADC (in the person of Mr. Frank Braddley) to ascertain the level 
of accuracy previously obtained by RADC in range-rate measurements. No 
precise information was obtained, though it was indicated that the accuracy 
was subject to question. For this reason, reference was again made to the 
raw data to determine whether variance data earlier assumed for range-rate data 
were reasonable. The results of this review indicated that these earlier 
data should be replaced with an estimate of variance on the order of .0001 
Km</sec®. Therefore, these modifications to the data enumerated in the 
previous sections have been made. 


At this point, refer to the sample problem presented on subsequent pages 
(pass 6069 of ECHO II) and note that 


1) All data utilized in the program are printed for reference during 
checkout. 
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2) The operation of the trajectory, integration and tracking groups 
is demonstrated. 


3) The behavior of the filter is demonstrated. 


Now, note that after the trace of the covariance matrix for estimation errors 
has been reduced (the initial errors are quite large since they have propagated 
for approximately 8 days) and orbital elements begin to be printed for the 
rectified (corrected) orbit (U.T. = 56616. sec),there is a steady change in the 
elements "a" and "e." The first printed values lie somewhat below those 
predicted by Goddard (about 1% in semi-major axis and 10% in eccentricity), 

but as subsequent points are processed the solution is attained, passed, 
attained again, and finally passed for the second time. 


As is apparent, the process is tending to converge, though the rate of 
convergence is less than originally anticipated. The observed rate of con- 
vergence is believed to be affected by two distinct factors. First, the 
nature of the doppler data itself will preclude exact convergence and may 
introduce an oscillatory motion in the solution. Second, the "memory" of 
the system is slight. In the way of explanation, the history (effect ) of 
all previous data points is carried in the covariance matrix for the 
estimation errors. This matrix in turn was allowed to "deteriorate" for almost 
eight days following the last "fix" on the satellite in the sample problem 
though, in reality, more recent estimates of this matrix probably existed. 
Thus, while the trajectory itself could be updated with reasonable accuracy 
from the prescribed epoch of 20 April 1965, it is unlikely that a "true" 
data reduction problem would be subjected to errors as large as those 
utilized in the sample except for the initial phases of the problem (i.e., 
immediately following injection). Rather, it is expected that the estimation 
errors at the epoch of 27 April 1965 were probably similar to those con- 
structed for the epoch of 20 April 1965. This change in the sample problem 
would reduce terms in the assumed covariance matrix by terms ranging from 
103 to 104, and would reduce the effect of the differences in the observed 
and computed values of the data by the same factors, thus allowing the 
trajectory of ECHO II to be corrected with precision. To illustrate this 
conclusion, the reduced data for anepoch near the final data point is 
presented in the following table for this variation of the sample problen. 


TABLE 2 


_FLLIPTIC @RBIT ELEMENTS 


_SEMIMAJOR AXIS 
ECCENTRICITY 

_TRUE ANGMOLY DATE (DEG) 
R ASC SF ASC NGONDF (DFG) 
ARG GF PERIAPSE (DEG) = 
GRBIT INCLINATIGN (NFG) 
OFL NGDE PER REV (RAD) 


DEL APSF PER REV (RAD) 


= AKM) O.T5307251F 04 
0. 24977477E-O1 
= __0214263121F OL 
0.24846007E O02 
_0+14930505E 02 

0.81473733E 02 
~0. 1OB8A1821E-02 


~0.32604282E-02 
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As is also apparent in the sample problem, the limit of the convergence 
process is not identical to the Goddard data. This fact is not alarming under 
the circumstances, due to the gross nature of the assumptions regarding the 
statistics of the problem, the poor quality of some of the data, and the 
small number of points processed. (The entire »ass was reduced to approxi- 
mately 45 points.) It is fully expected that a more complete sample problem 
involving more precisely reduced raw data will converge to the desired 
solution. 
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5.0 CONCLUSIONS AND RECOMMENDATIONS 


The program described in the previous pages is an accurate and efficient 
tool for estimating the trajectory of a satellite traveling about the earth. 
There are, however, several features which should be incorporated in the 
program to extend the number of applications and to further enhance the 
accuracy. The most important of these are itemized below (not in the order 
of importance). 


1) The present program prints at intervals determined by the step 
size employed for the numerical integration. This mode of operation 
is completely adequate for applications of the type of the sample 
problem. However, should the program be applied to general trajectory 
problems (a simple task implemented by replacing the present sub- 
routine DATAPE with a dummy which reads an end time and the sub- 
rontine FILTER with a routine which terminates oneration if called), 
prints at regular prescribed intervals would normally be required, 
thus, an alternate logic for this application is recommended. 


2) ‘the present logic employed in the integration cycle involves 
restarting the integration with the most recent values of r and ¥ 
at those times when the step size has been changed. An imnroved 
logic would be to use the undated differences and sums for an 
epoch 3 steps before the most recent point to reevaluate the fF and 
¥ at the time in the past using the central difference formula for 
integration. This procedure would assure that errors introduced 
in the extrapolation of ¥ and ¥ could be controlled to a higher 
degree. 


3) The present starter routine for the numerical integration is a 4th 
order Runge-Kutta. This routine provides a series of values of fF, 
¥ and @ which are differenced and summed to construct the table 
required ‘in the Gauss-Jackson stepping. However, the Gauss routine 
operating with the predictor cycle done is a more precise routine 
than the R-K starter, Thus, if extreme precision is desired, the 
difference table provided by the start cycle must be corrected. 
This correction could be achieved by assuming the 6th difference 
is constant, constructing the missing differences back to To; Vos 
and employing the central difference formula for reestimating the : 
values of ¥ and V to be used for the difference table. ~ 


1) The constants of the math model (the coefficients of the potertial 
function, the gravitational constant, the quantity Cpa/M, etc.) 
could be estimated at the same time the state is being differentially 
corrected to adjust for slight bias errors. This mode of operation 
would assure that the trajectory could be estimated to a higher 
degree of precision than presently obtainable. Alternate logic 
would, of course, be required to determine when sufficient information 
w2s avilable in the data to rake such a process convergent. 
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The present Josic makes no at‘ empt to correct the raw data for any 
known biases, for biases of the trne displaved in Figures 2 and 5 
(discussed ir the sample problom), for time of sicnal propagation, 

for refraction, ete. Thus, these effects must be reflected in the 

data prior to processing in the PREPROCESSOR. Slight modification 

of the preliminary logic could, however, effect this correction and 

save considerable effort at the recording station. Thus, if corrections 
to the data are worth; of inclusion, these changes must be made. 


6) Some editing of the raw data is prescntly accomplished (the observa- 
tiors are made "continuous" functions of time; o.g., if the data 
are increasing toward 27 and a value of slightly more than zero is 
encountered, adjustment will bc made to assure that all data fall 
near some contirucus curve before smoothing). However, no editing 
to determine if time is recorded properly is attempted. Several 
tests of this latter nature should be devised to nmiscard the bad 
data points so that the possibility of effectine the smoothed 
data files can te disregarded (sce the discussion of the preliminary 
processor in the sample problem). 


7) The present system of providing the input data misht conceivably 
be revised to simplify operation. Particular attention should be 
directed to the manner in which time is input. This data could 
he provided in a more straightforward manner in the form of vear, 


month, day, hour, second, or ‘lian date. 


Other modifications would of course be required to allow lunar and planetary 
trajectories to be analyzed. These modifications would all stem from the 
inclusior of an accurate ephemeris routine and the logic necessary to produce 
information in other coordinate systems. Should this task be adjudged worthy 
of attention, simultaneous effort should be directed toward the inclusion of 
midcourse guidance logic so that simulations of true trajectories and studies 
of trajectory shapins can be effected with efficiency. 


The application of this program to real problems (such as that which 
was attempted in the sample) requires precise data for the station, for the 
station errors, for the noise, and for the observations (this latter point. 
was discussed in items 5 and 6). However, no means is provided within the 
logic as a check on these data. Thus, it is considered essential that the 
Vacility be periodically calibrated against the remainder of the tracking 
network, and that the results of this calibration he whilized in preference 
to assumed data of the types built inte subroutine TYPUT. 


Firelly, i: is recomnended thet enother sample problem be prepared. 
this new sample should utilize precisely reduced observed data (of all six 
types permissible) from several tracking stations (acquired over a period of 
several orbits), precise station data, precise noise data, and precise error 
data for the station's position. Only in this manner can a final check be 
performed. 
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APPENDIX I 
Coenen 


CANMgn data storase is affected through the facility of labeled (non- 
bla’ k) and unlabeled (blank) ogiMen , A detailed descri;ytion of the contents 
of both storage facilities is contained within this section. 


Biank C@MK@N Map. 


A blank (unlabeled) C¢MMEN section is utilized as the principal means 
of communication between the various subprorsrams. This section, specified 
by the 525 element DATA array, is composed of the following five subarrays. 


CON: A 15 element subarray, loaded at DATA (1), containing program 
constants and conversion factors. 


SAT: A 20 element subarray, loaded at DATA (16), containing data 
pertaining to the satellite at the initial epoch. 


SDA; A 250 element subarray, loaded at DATA (36), containing data 
pertaining to the tracking stations. 


STT. A 105 element subarray, loaded at DATA (286), containing 
statistical information and releted data for the satellite. 


WRK: A 135 element subarray loaded at DATA (391). The primary 
function of the array is to act asa "scratch rarer" or 
working array for communication between the various subroutines, 


The following map, which describes the bl ink CAMNAN region in detail, 
utilizes standard FORTRAN nomenclature with a single exception. When a 
description refers to a series of locations in a DATA subarray (usually 
time), an abbreviated notati-n is used, in that, reference to the two cells 
WRK (50) and WRK (51) is written WRK (50,51). From the context, it will 
be apparent that WRK (50,51) is not a single element of a double sub- 
scripted array. 
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Labeled CYMIM~N MAP 


In audition to blank C@MM@N, COMNG@N storage utilized by the pro--rran 
includes three regions of labeled C@MiGN. These latter regions, containing 
tabulated atmospheric and ephemeris data, are set at load time by means of 
a BLOCK DATA subprogram. These regions are explicitly accessed by sub- 
routines ATMS and EPHEM exclusively. A description of the data stored in 
the labeled COMMeN regions is tabulated below. 
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APPENDIX II 


{ PROGRAM INPUT/OUTPUT TAPES 


The program logic has been constructed in such a manner that variable 
input/output tape numbers can be employed so that operation on any system with 
a basic FORTRAN capability will be possible with a minimum of effort, This 
mode of operation has been provided by placing the desired tape numbers in 
COMMON [CON(16), CON(17)] while loading subroutine INPUT and calling for them 
in each routine requiring input or output. 





a 
Unfortunately, due to problems of checkout with the present NAASYS system, 
it was necessary to revise this logic and insert the actual tape units utilized 
in each routine. 
INPUT tape unit 5 
OUT PUT tape unit 6 
SPECLAL tape unit 9 
The linkage to COMMON was not, nowever, removed. Thus, if suffziciant CORE 
exists to allow communication with each tape unit in this mode, the capability 
should be reinstated to assure complete compatibility of the program with all 
; FORTRAN systems and to avoid problems associated with updating or changing the 
system, 
h For the sake of reference, the locations of all input/output statements 
in the program are listed helow, 
SUBROUTINE 
an 
MAIN 
INPUT 3910, 3931 
REED 0410 
se SICRD 0058 
DADUMP 0250, 0330, 0370, 0450, 0490, 0500, 0580, 0620, 0700, 0740, 0795 
ATMS 
TRAK 1022, 1480, 1539, 1550 
UPSTAT 0370, 0400, 0420, 0660, 0770 
DATAPE 0720, 0880 
ATANS 
oa 
4 TT-1 





nreme emanemrans, 


\~ 


The preceding comments pertain only to those systems for which the use of 
{  vontrol cards (in the binary deck) to equivalence the addresses of the tape 

i units cequested and those actually utilized is not permitted. However, if 

{ su.h cperation is possible, this "fix" is by far the simplest remedy. 
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